forked from janelia-flyem/dvid
/
geometry.go
651 lines (577 loc) · 17 KB
/
geometry.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
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
package dvid
import (
"fmt"
"strings"
"sync"
)
type Bounder interface {
// StartPoint returns the offset to first point of data.
StartPoint() Point
// EndPoint returns the last point.
EndPoint() Point
}
// Geometry describes the shape, size, and position of data in the DVID volume.
type Geometry interface {
Bounder
// DataShape describes the shape of the data.
DataShape() DataShape
// Size returns the extent in each dimension.
Size() Point
// NumVoxels returns the number of voxels within this space.
NumVoxels() int64
String() string
}
// Extents holds the extents of a volume in both absolute voxel coordinates
// and lexicographically sorted chunk indices.
type Extents struct {
MinPoint Point
MaxPoint Point
MinIndex ChunkIndexer
MaxIndex ChunkIndexer
pointMu sync.Mutex
indexMu sync.Mutex
}
// Duplicate creates a duplicate Extents.
func (ext *Extents) Duplicate() Extents {
var dup Extents
dup.MinPoint = ext.MinPoint.Duplicate()
dup.MaxPoint = ext.MaxPoint.Duplicate()
dup.MinIndex = ext.MinIndex.DuplicateChunkIndexer()
dup.MaxIndex = ext.MaxIndex.DuplicateChunkIndexer()
return dup
}
// --- dvid.Bounder interface ----
// StartPoint returns the offset to first point of data.
func (ext *Extents) StartPoint() Point {
return ext.MinPoint
}
// EndPoint returns the last point.
func (ext *Extents) EndPoint() Point {
return ext.MaxPoint
}
// --------
// AdjustPoints modifies extents based on new voxel coordinates in concurrency-safe manner.
func (ext *Extents) AdjustPoints(pointBeg, pointEnd Point) bool {
ext.pointMu.Lock()
defer ext.pointMu.Unlock()
var changed bool
if ext.MinPoint == nil {
ext.MinPoint = pointBeg
changed = true
} else {
ext.MinPoint, changed = ext.MinPoint.Min(pointBeg)
}
if ext.MaxPoint == nil {
ext.MaxPoint = pointEnd
changed = true
} else {
ext.MaxPoint, changed = ext.MaxPoint.Max(pointEnd)
}
return changed
}
// AdjustIndices modifies extents based on new block indices in concurrency-safe manner.
func (ext *Extents) AdjustIndices(indexBeg, indexEnd ChunkIndexer) bool {
ext.indexMu.Lock()
defer ext.indexMu.Unlock()
var changed bool
if ext.MinIndex == nil {
ext.MinIndex = indexBeg
changed = true
} else {
ext.MinIndex, changed = ext.MinIndex.Min(indexBeg)
}
if ext.MaxIndex == nil {
ext.MaxIndex = indexEnd
changed = true
} else {
ext.MaxIndex, changed = ext.MaxIndex.Max(indexEnd)
}
return changed
}
// GetNumBlocks returns the number of n-d blocks necessary to cover the given geometry.
func GetNumBlocks(geom Geometry, blockSize Point) int {
startPt := geom.StartPoint()
size := geom.Size()
numBlocks := 1
for dim := uint8(0); dim < geom.Size().NumDims(); dim++ {
blockLength := blockSize.Value(dim)
startMod := startPt.Value(dim) % blockLength
length := size.Value(dim) + startMod
blocks := length / blockLength
if length%blockLength != 0 {
blocks++
}
numBlocks *= int(blocks)
}
return numBlocks
}
type Dimension struct {
Name string
Units string
beg, end int32
}
var (
// XY describes a 2d rectangle of voxels that share a z-coord.
XY = DataShape{3, []uint8{0, 1}}
// XZ describes a 2d rectangle of voxels that share a y-coord.
XZ = DataShape{3, []uint8{0, 2}}
// YZ describes a 2d rectangle of voxels that share a x-coord.
YZ = DataShape{3, []uint8{1, 2}}
// Arb describes a 2d rectangle of voxels with arbitrary 3d orientation.
Arb = DataShape{3, nil}
// Vol3d describes a 3d volume of voxels
Vol3d = DataShape{3, []uint8{0, 1, 2}}
)
// DataShape describes the number of dimensions and the ordering of the dimensions.
type DataShape struct {
dims uint8
shape []uint8
}
const DataShapeBytes = 7
// BytesToDataShape recovers a DataShape from a series of bytes.
func BytesToDataShape(b []byte) (s DataShape, err error) {
if b == nil {
err = fmt.Errorf("Cannot convert nil to DataShape!")
return
}
if len(b) != DataShapeBytes {
err = fmt.Errorf("Cannot convert %d bytes to DataShape", len(b))
return
}
s = DataShape{dims: uint8(b[0])}
shapeLen := int(b[1])
s.shape = make([]uint8, shapeLen)
for i := 0; i < shapeLen; i++ {
s.shape[i] = b[i+2]
}
return
}
// AxisName returns common axis descriptions like X, Y, and Z for a shapes dimensions.
func (s DataShape) AxisName(axis uint8) string {
if int(axis) >= len(s.shape) {
return "Unknown"
}
switch s.shape[axis] {
case 0:
return "X"
case 1:
return "Y"
case 2:
return "Z"
default:
return fmt.Sprintf("Dim %d", axis)
}
}
// Bytes returns a fixed length byte representation that can be used for keys.
// Up to 5-d shapes can be used.
func (s DataShape) Bytes() []byte {
b := make([]byte, DataShapeBytes)
b[0] = byte(s.dims)
b[1] = byte(len(s.shape))
for i := 0; i < len(s.shape); i++ {
b[i+2] = s.shape[i]
}
return b
}
// TotalDimensions returns the full dimensionality of space within which there is this DataShape.
func (s DataShape) TotalDimensions() int8 {
return int8(s.dims)
}
// ShapeDimensions returns the number of dimensions for this shape.
func (s DataShape) ShapeDimensions() int8 {
if s.shape == nil {
return 0
}
return int8(len(s.shape))
}
// ShapeDimension returns the axis number for a shape dimension.
func (s DataShape) ShapeDimension(axis uint8) (uint8, error) {
if s.shape == nil {
return 0, fmt.Errorf("Cannot request ShapeDimension from nil DataShape")
}
if len(s.shape) <= int(axis) {
return 0, fmt.Errorf("Illegal dimension requested from DataShape: %d", axis)
}
return s.shape[axis], nil
}
// GetSize2D returns the width and height of a 2D shape given a n-D size.
func (s DataShape) GetSize2D(size SimplePoint) (width, height int32, err error) {
if len(s.shape) != 2 {
err = fmt.Errorf("Can't get 2D size from a non-2D shape: %s", s)
return
}
width = size.Value(s.shape[0])
height = size.Value(s.shape[1])
return
}
// GetFloat2D returns elements of a N-d float array given a 2d shape.
func (s DataShape) GetFloat2D(fslice NdFloat32) (x, y float32, err error) {
if len(s.shape) != 2 {
err = fmt.Errorf("Can't get 2D data from a non-2D shape: %s", s)
return
}
x = fslice[s.shape[0]]
y = fslice[s.shape[1]]
return
}
// ChunkPoint3d returns a chunk point where the XY is determined by the
// type of slice orientation of the DataShape, and the Z is the non-chunked
// coordinate. This is useful for tile generation where you have 2d tiles
// in a 3d space.
func (s DataShape) ChunkPoint3d(p, size Point) (ChunkPoint3d, error) {
if len(s.shape) != 2 {
return ChunkPoint3d{}, fmt.Errorf("Can't get process slice from a non-2D shape: %s", s)
}
if s.dims != 3 {
return ChunkPoint3d{}, fmt.Errorf("ChunkPoint3d() can only be called on 3d points!")
}
chunkable, ok := p.(Chunkable)
if !ok {
return ChunkPoint3d{}, fmt.Errorf("ChunkPoint3d() requires Chunkable point.")
}
chunk := chunkable.Chunk(size)
switch {
case s.Equals(XY):
return ChunkPoint3d{chunk.Value(0), chunk.Value(1), p.Value(2)}, nil
case s.Equals(XZ):
return ChunkPoint3d{chunk.Value(0), chunk.Value(2), p.Value(1)}, nil
case s.Equals(YZ):
return ChunkPoint3d{chunk.Value(1), chunk.Value(2), p.Value(0)}, nil
default:
return ChunkPoint3d{}, fmt.Errorf("ChunkPoint3d() can only be run on slices: given %s", s)
}
}
// PlaneToChunkPoint3d returns a chunk point corresponding to the given point on the DataShape's
// plane. If DataShape is not a plane, returns an error.
func (s DataShape) PlaneToChunkPoint3d(x, y int32, offset, size Point) (ChunkPoint3d, error) {
if len(s.shape) != 2 {
return ChunkPoint3d{}, fmt.Errorf("Can't get plane point from a non-2D shape: %s", s)
}
if s.dims != 3 {
return ChunkPoint3d{}, fmt.Errorf("PlaneToChunkPoint3d() requires 3D shape: %s", s)
}
var p Point3d
switch {
case s.Equals(XY):
p = Point3d{x + offset.Value(0), y + offset.Value(1), offset.Value(2)}
chunkPt := p.Chunk(size)
return ChunkPoint3d{chunkPt.Value(0), chunkPt.Value(1), p[2]}, nil
case s.Equals(XZ):
p = Point3d{x + offset.Value(0), offset.Value(1), y + offset.Value(2)}
chunkPt := p.Chunk(size)
return ChunkPoint3d{chunkPt.Value(0), p[1], chunkPt.Value(2)}, nil
case s.Equals(YZ):
p = Point3d{offset.Value(0), x + offset.Value(1), y + offset.Value(2)}
chunkPt := p.Chunk(size)
return ChunkPoint3d{p[0], chunkPt.Value(1), chunkPt.Value(2)}, nil
default:
return ChunkPoint3d{}, fmt.Errorf("ChunkPoint3d() can only be run on slices: given %s", s)
}
}
// Duplicate returns a duplicate of the DataShape.
func (s DataShape) Duplicate() DataShape {
return DataShape{dims: s.dims, shape: append([]uint8{}, s.shape...)}
}
// Equals returns true if the passed DataShape is identical.
func (s DataShape) Equals(s2 DataShape) bool {
if s.dims == s2.dims && len(s.shape) == len(s2.shape) {
for i, dim := range s.shape {
if s2.shape[i] != dim {
return false
}
}
return true
}
return false
}
func (s DataShape) String() string {
switch {
case s.Equals(XY):
return "XY slice"
case s.Equals(XZ):
return "XZ slice"
case s.Equals(YZ):
return "YZ slice"
case s.Equals(Arb):
return "slice with arbitrary orientation"
case s.Equals(Vol3d):
return "3d volume"
case s.dims > 3:
return "n-D volume"
default:
return "Unknown shape"
}
}
// DataShapes are a slice of DataShape.
type DataShapes []DataShape
// GetShapes returns DataShapes from a string where each shape specification
// is delimited by a separator. If the key is not found, nil is returned.
func (c Config) GetShapes(key, separator string) ([]DataShape, error) {
s, found, err := c.GetString(key)
if err != nil {
return nil, err
}
if !found {
return nil, nil
}
parts := strings.Split(s, separator)
shapes := []DataShape{}
for _, part := range parts {
shape, err := DataShapeString(part).DataShape()
if err != nil {
return nil, err
}
shapes = append(shapes, shape)
}
return shapes, nil
}
// String for specifying a slice orientation or subvolume
type DataShapeString string
// List of strings associated with shapes up to 3d
var dataShapeStrings = map[string]DataShape{
"xy": XY,
"xz": XZ,
"yz": YZ,
"vol": Vol3d,
"arb": Arb,
"0_1": XY,
"0_2": XZ,
"1_2": YZ,
"0_1_2": Vol3d,
"0,1": XY,
"0,2": XZ,
"1,2": YZ,
"0,1,2": Vol3d,
}
// ListDataShapes returns a slice of shape names
func ListDataShapes() (shapes []string) {
shapes = []string{}
for key, _ := range dataShapeStrings {
shapes = append(shapes, string(key))
}
return
}
// DataShape returns the data shape constant associated with the string.
func (s DataShapeString) DataShape() (shape DataShape, err error) {
shape, found := dataShapeStrings[strings.ToLower(string(s))]
if !found {
err = fmt.Errorf("Unknown data shape specification (%s)", s)
}
return
}
// Returns the image size necessary to compute an isotropic slice of the given dimensions.
// If isotropic is false, simply returns the original slice geometry. If isotropic is true,
// uses the higher resolution dimension.
func Isotropy2D(voxelSize NdFloat32, geom Geometry, isotropic bool) (Geometry, error) {
if !isotropic {
return geom, nil
}
// Get the voxel resolutions for this particular slice orientation
resX, resY, err := geom.DataShape().GetFloat2D(voxelSize)
if err != nil {
return nil, err
}
if resX == resY {
return geom, nil
}
srcW := geom.Size().Value(0)
srcH := geom.Size().Value(1)
var dstW, dstH int32
if resX < resY {
// Use x resolution for all pixels.
dstW = srcW
dstH = int32(float32(srcH)*resX/resY + 0.5)
} else {
dstH = srcH
dstW = int32(float32(srcW)*resY/resX + 0.5)
}
// Make altered geometry
slice, ok := geom.(*OrthogSlice)
if !ok {
return nil, fmt.Errorf("can only handle isotropy for orthogonal 2d slices")
}
dstSlice := slice.Duplicate()
dstSlice.SetSize(Point2d{dstW, dstH})
return dstSlice, nil
}
// ---- Geometry implementations ------
// Subvolume describes a 3d box Geometry. The "Sub" prefix emphasizes that the
// data is usually a smaller portion of the volume held by the DVID datastore.
// Note that the 3d coordinate system is assumed to be a Right-Hand system like OpenGL.
type Subvolume struct {
shape DataShape
offset Point
size Point
}
// NewSubvolumeFromStrings returns a Subvolume given string representations of
// offset ("0,10,20") and size ("250,250,250").
func NewSubvolumeFromStrings(offsetStr, sizeStr, sep string) (*Subvolume, error) {
offset, err := StringToPoint(offsetStr, sep)
if err != nil {
return nil, err
}
if offset.NumDims() != 3 {
return nil, fmt.Errorf("Offset must be 3d coordinate, not %d-d coordinate", offset.NumDims())
}
size, err := StringToPoint(sizeStr, sep)
if err != nil {
return nil, err
}
if size.NumDims() != 3 {
return nil, fmt.Errorf("Size must be 3 (not %d) dimensions", size.NumDims())
}
return NewSubvolume(offset, size), nil
}
// NewSubvolume returns a Subvolume given a subvolume's origin and size.
func NewSubvolume(offset, size Point) *Subvolume {
return &Subvolume{Vol3d, offset, size}
}
func (s *Subvolume) DataShape() DataShape {
return Vol3d
}
func (s *Subvolume) Size() Point {
return s.size
}
func (s *Subvolume) NumVoxels() int64 {
if s == nil || s.size.NumDims() == 0 {
return 0
}
voxels := int64(s.size.Value(0))
for dim := uint8(1); dim < s.size.NumDims(); dim++ {
voxels *= int64(s.size.Value(dim))
}
return voxels
}
func (s *Subvolume) StartPoint() Point {
return s.offset
}
func (s *Subvolume) EndPoint() Point {
return s.offset.Add(s.size.Sub(Point3d{1, 1, 1}))
}
func (s *Subvolume) String() string {
return fmt.Sprintf("%s %s at offset %s", s.shape, s.size, s.offset)
}
// NewIndexZYXIterator returns an iterator that can move across subvolume.
func (s *Subvolume) NewIndexZYXIterator(chunkSize Point) (IndexIterator, error) {
// Setup traversal
begVoxel, ok := s.StartPoint().(Chunkable)
if !ok {
return nil, fmt.Errorf("Subvolume %s StartPoint() cannot handle Chunkable points.", s)
}
endVoxel, ok := s.EndPoint().(Chunkable)
if !ok {
return nil, fmt.Errorf("Subvolume %s EndPoint() cannot handle Chunkable points.", s)
}
begBlock, ok := begVoxel.Chunk(chunkSize).(ChunkPoint3d)
if !ok {
return nil, fmt.Errorf("Subvolume %s StartPoint() is not a 3d chunk", s)
}
endBlock, ok := endVoxel.Chunk(chunkSize).(ChunkPoint3d)
if !ok {
return nil, fmt.Errorf("Subvolume %s EndPoint() is not a 3d chunk", s)
}
return NewIndexZYXIterator(begBlock, endBlock), nil
}
// OrthogSlice is a 2d rectangle orthogonal to two axis of the space that is slices.
// It fulfills a Geometry interface.
type OrthogSlice struct {
shape DataShape
offset Point
size Point2d
endPoint Point
}
// NewSliceFromStrings returns a Geometry object for a XY, XZ, or YZ slice given
// a data shape string, offset ("0,10,20"), and size ("250,250").
func NewSliceFromStrings(str DataShapeString, offsetStr, sizeStr, sep string) (Geometry, error) {
shape, err := str.DataShape()
if err != nil {
return nil, err
}
offset, err := StringToPoint(offsetStr, sep)
if err != nil {
return nil, err
}
if offset.NumDims() != 3 {
return nil, fmt.Errorf("Offset must be 3d coordinate, not %d-d coordinate", offset.NumDims())
}
// Enforce that size string is 2d since this is supposed to be a slice.
ndstring, err := StringToNdString(sizeStr, sep)
if err != nil {
return nil, err
}
size, err := ndstring.Point2d()
if err != nil {
return nil, err
}
return NewOrthogSlice(shape, offset, size)
}
// NewOrthogSlice returns an OrthogSlice of chosen orientation, offset, and size.
func NewOrthogSlice(s DataShape, offset Point, size Point2d) (Geometry, error) {
if offset.NumDims() != s.dims {
return nil, fmt.Errorf("NewOrthogSlice: offset dimensionality %d != shape %d",
offset.NumDims(), s.dims)
}
if s.shape == nil || len(s.shape) != 2 {
return nil, fmt.Errorf("NewOrthogSlice: shape not properly specified")
}
xDim := s.shape[0]
if xDim >= s.dims {
return nil, fmt.Errorf("NewOrthogSlice: X dimension of slice (%d) > # avail dims (%d)",
xDim, s.dims)
}
yDim := s.shape[1]
if yDim >= s.dims {
return nil, fmt.Errorf("NewOrthogSlice: Y dimension of slice (%d) > # avail dims (%d)",
yDim, s.dims)
}
settings := map[uint8]int32{
xDim: offset.Value(xDim) + size[0] - 1,
yDim: offset.Value(yDim) + size[1] - 1,
}
geom := &OrthogSlice{
shape: s,
offset: offset.Duplicate(),
size: size,
endPoint: offset.Modify(settings),
}
return geom, nil
}
func (s *OrthogSlice) SetSize(size Point2d) {
s.size = size
xDim := s.shape.shape[0]
yDim := s.shape.shape[1]
settings := map[uint8]int32{
xDim: s.offset.Value(xDim) + size[0] - 1,
yDim: s.offset.Value(yDim) + size[1] - 1,
}
s.endPoint = s.offset.Modify(settings)
}
func (s OrthogSlice) Duplicate() OrthogSlice {
return OrthogSlice{
shape: s.shape.Duplicate(),
offset: s.offset.Duplicate(),
size: s.size.Duplicate().(Point2d),
endPoint: s.endPoint.Duplicate(),
}
}
// --- Geometry interface -----------
func (s OrthogSlice) DataShape() DataShape {
return s.shape
}
func (s OrthogSlice) Size() Point {
return s.size
}
func (s OrthogSlice) NumVoxels() int64 {
x := int64(s.size[0])
y := int64(s.size[1])
return x * y
}
func (s OrthogSlice) StartPoint() Point {
return s.offset
}
func (s OrthogSlice) EndPoint() Point {
return s.endPoint
}
func (s OrthogSlice) String() string {
return fmt.Sprintf("%s @ offset %s, size %s", s.shape, s.offset, s.size)
}