github.com/janelia-flyem/dvid@v1.0.0/datatype/common/labels/compressed_old.go (about)

     1  // This file implements the initial version of Google Block compression by Steve Plaza.
     2  // It is retained for benchmarks as a proxy for standard Neuroglancer-based Block compression.
     3  
     4  package labels
     5  
     6  import (
     7  	"fmt"
     8  	"math"
     9  	"runtime"
    10  
    11  	"github.com/janelia-flyem/dvid/dvid"
    12  )
    13  
    14  // oldCompress takes an input 3D label array of a given volume dimensions and compresses
    15  // it using the segmentation format defined in Neuroglancer.  VolSize must be tileable by
    16  // 8x8x8 block size.
    17  func oldCompress(b []byte, volsize dvid.Point) ([]byte, error) {
    18  	uintSlice, err := dvid.AliasByteToUint64(b)
    19  	if err != nil {
    20  		return nil, err
    21  	}
    22  	seg, err := oldCompressUint64(uintSlice, volsize)
    23  	if err != nil {
    24  		return nil, err
    25  	}
    26  	runtime.KeepAlive(&b)
    27  	return dvid.AliasUint32ToByte(seg), nil
    28  }
    29  
    30  // oldCompressUint64 takes an input 3D label array of a given volume dimensions and compresses
    31  // it using the segmentation format defined in Neuroglancer.  VolSize must be tileable by
    32  // 8x8x8 block size.  The output data is an array of 32 bit numbers. TODO: reuse table
    33  // encoding to improve compression.
    34  func oldCompressUint64(input []uint64, volsize dvid.Point) ([]uint32, error) {
    35  	BLKSIZE := uint(8)
    36  
    37  	xsize := uint(volsize.Value(0))
    38  	ysize := uint(volsize.Value(1))
    39  	zsize := uint(volsize.Value(2))
    40  	gx := xsize / BLKSIZE
    41  	gy := ysize / BLKSIZE
    42  	gz := zsize / BLKSIZE
    43  	if xsize%BLKSIZE > 0 || ysize%BLKSIZE > 0 || zsize%BLKSIZE > 0 {
    44  		return nil, fmt.Errorf("volume must be a multiple of the block size")
    45  	}
    46  
    47  	// must add initial 4 byte to designate as a header
    48  	// for the compressed data for neuroglancer
    49  
    50  	// 64 bit headers for each 8x8x8 block and pre-allocate some data based on expected data size
    51  	globaloffset := 0 // !! set to 1 if header used by neuroglancer which needs to know first byte of multi-channel compression
    52  
    53  	datagoogle := make([]uint32, gx*gy*gz*2+uint(globaloffset), xsize*ysize*zsize*2/10)
    54  	//datagoogle[0] = byte(1) // !! compressed data starts after first 4 bytes for neuroglancer
    55  
    56  	// everything is written out little-endian
    57  	for gziter := uint(0); gziter < gz; gziter++ {
    58  		for gyiter := uint(0); gyiter < gy; gyiter++ {
    59  			for gxiter := uint(0); gxiter < gx; gxiter++ {
    60  				unique_vals := make(map[uint64]uint32)
    61  				unique_list := make([]uint64, 0)
    62  
    63  				currpos := (gziter*BLKSIZE*(xsize*ysize) + gyiter*BLKSIZE*xsize + gxiter*BLKSIZE)
    64  
    65  				// extract unique values in the 8x8x8 block
    66  				for z := uint(0); z < BLKSIZE; z++ {
    67  					for y := uint(0); y < BLKSIZE; y++ {
    68  						for x := uint(0); x < BLKSIZE; x++ {
    69  							if _, ok := unique_vals[input[currpos]]; !ok {
    70  								unique_vals[input[currpos]] = 0
    71  								unique_list = append(unique_list, input[currpos])
    72  							}
    73  							currpos += 1
    74  						}
    75  						currpos += (xsize - BLKSIZE)
    76  					}
    77  					currpos += (xsize*ysize - (xsize * (BLKSIZE)))
    78  				}
    79  				// write out mapping
    80  				for pos, val := range unique_list {
    81  					unique_vals[val] = uint32(pos)
    82  				}
    83  
    84  				// write-out compressed data
    85  				encodedBits := uint32(math.Ceil(math.Log2(float64(len(unique_vals)))))
    86  				switch {
    87  				case encodedBits == 0, encodedBits == 1, encodedBits == 2:
    88  				case encodedBits <= 4:
    89  					encodedBits = 4
    90  				case encodedBits <= 8:
    91  					encodedBits = 8
    92  				case encodedBits <= 16:
    93  					encodedBits = 16
    94  				}
    95  
    96  				// starting location for writing out data
    97  				compressstart := len(datagoogle) // in 4-byte units
    98  				// number of bytes to add (encode bytes + table size of 8 byte numbers)
    99  				addedInts := uint32(encodedBits*uint32(BLKSIZE*BLKSIZE*BLKSIZE)/8)/4 + uint32(len(unique_vals)*2) // will always be a multiple of 4 bytes
   100  				bitspot := uint(len(datagoogle) * 32)
   101  				datagoogle = append(datagoogle, make([]uint32, addedInts)...)
   102  
   103  				// do not need to write-out anything if there is only one entry
   104  				if encodedBits > 0 {
   105  					currpos := (gziter*BLKSIZE*(xsize*ysize) + gyiter*BLKSIZE*xsize + gxiter*BLKSIZE)
   106  					for z := uint32(0); z < uint32(BLKSIZE); z++ {
   107  						for y := uint32(0); y < uint32(BLKSIZE); y++ {
   108  							for x := uint32(0); x < uint32(BLKSIZE); x++ {
   109  								mappedval := unique_vals[input[currpos]]
   110  								bitshift := bitspot % 32
   111  								bytespot := bitspot / 32
   112  								datagoogle[bytespot] |= (mappedval << bitshift)
   113  
   114  								bitspot += uint(encodedBits)
   115  								currpos += 1
   116  							}
   117  							currpos += (xsize - BLKSIZE)
   118  						}
   119  						currpos += (xsize*ysize - (xsize * (BLKSIZE)))
   120  					}
   121  				}
   122  
   123  				// write-out lookup table
   124  				tablestart := bitspot / 32 // in 4-byte units
   125  				for _, val := range unique_list {
   126  					datagoogle[bitspot/32] = uint32(val)
   127  					bitspot += 32
   128  					datagoogle[bitspot/32] = uint32(val >> 32)
   129  					bitspot += 32
   130  				}
   131  
   132  				// write-out block header
   133  				// 8 bytes per header entry
   134  				headerpos := (gziter*(gy*gx)+gyiter*gx+gxiter)*2 + uint(globaloffset) // shift start by global offset
   135  
   136  				// write out lookup table start
   137  				tablestart -= uint(globaloffset) // relative to the start of the compressed data
   138  				datagoogle[headerpos] = uint32(tablestart)
   139  
   140  				// write out number of encoded bits
   141  				datagoogle[headerpos] |= (encodedBits << 24)
   142  				headerpos++
   143  
   144  				// write out block compress start
   145  				compressstart -= (globaloffset) // relative to the start of the compressed data
   146  				datagoogle[headerpos] = uint32(compressstart)
   147  				headerpos++
   148  			}
   149  		}
   150  	}
   151  
   152  	return datagoogle, nil
   153  
   154  }
   155  
   156  // oldDecompress takes an input compressed array of 64-bit labels of a given volume dimensions
   157  // and decompresses it. VolSize must be tileable by 8x8x8 block size.  The output
   158  // data is a slice of bytes that represents packed uint64.
   159  func oldDecompress(b []byte, volsize dvid.Point) ([]byte, error) {
   160  	seg, err := dvid.AliasByteToUint32(b)
   161  	if err != nil {
   162  		return nil, err
   163  	}
   164  	labels, err := oldDecompressUint64(seg, volsize)
   165  	if err != nil {
   166  		return nil, err
   167  	}
   168  	return dvid.AliasUint64ToByte(labels), nil
   169  }
   170  
   171  // oldDecompressUint64 takes an input compressed array of 64-bit labels of a given volume dimensions
   172  // and decompresses it. VolSize must be tileable by 8x8x8 block size.  The output
   173  // data is an array of 64 bit numbers.
   174  func oldDecompressUint64(input []uint32, volsize dvid.Point) ([]uint64, error) {
   175  	BLKSIZE := uint(8)
   176  
   177  	xsize := uint(volsize.Value(0))
   178  	ysize := uint(volsize.Value(1))
   179  	zsize := uint(volsize.Value(2))
   180  	if xsize%BLKSIZE > 0 || ysize%BLKSIZE > 0 || zsize%BLKSIZE > 0 {
   181  		return nil, fmt.Errorf("volume must be a multiple of the block size")
   182  	}
   183  
   184  	// create output buffer
   185  	output := make([]uint64, xsize*ysize*zsize)
   186  
   187  	grid_size := [3]uint{xsize / BLKSIZE, ysize / BLKSIZE, zsize / BLKSIZE}
   188  	block := [3]uint{0, 0, 0}
   189  
   190  	for block[2] = 0; block[2] < grid_size[2]; block[2]++ {
   191  		for block[1] = 0; block[1] < grid_size[1]; block[1]++ {
   192  			for block[0] = 0; block[0] < grid_size[0]; block[0]++ {
   193  
   194  				block_offset := block[0] + grid_size[0]*(block[1]+grid_size[1]*block[2])
   195  
   196  				tableoffset := input[block_offset*2] & 0xffffff
   197  				encoded_bits := (input[block_offset*2] >> 24) & 0xff
   198  				encoded_value_start := input[block_offset*2+1]
   199  
   200  				// find absolute positions in output array
   201  				xmin := block[0] * BLKSIZE
   202  				xmax := xmin + BLKSIZE
   203  
   204  				ymin := block[1] * BLKSIZE
   205  				ymax := ymin + BLKSIZE
   206  
   207  				zmin := block[2] * BLKSIZE
   208  				zmax := zmin + BLKSIZE
   209  
   210  				bitmask := (1 << encoded_bits) - 1
   211  
   212  				for z := zmin; z < zmax; z++ {
   213  					for y := ymin; y < ymax; y++ {
   214  
   215  						outindex := (z*ysize+y)*xsize + xmin
   216  						bitpos := BLKSIZE * ((z-zmin)*(BLKSIZE) + (y - ymin)) * uint(encoded_bits)
   217  
   218  						for x := xmin; x < xmax; x++ {
   219  							bitshift := bitpos % 32
   220  
   221  							arraypos := bitpos / (32)
   222  							bitval := uint(0)
   223  							if encoded_bits > 0 {
   224  								bitval = (uint(input[uint(encoded_value_start)+arraypos]) >> bitshift) & uint(bitmask)
   225  							}
   226  							val := uint64(input[uint(tableoffset)+bitval*2])
   227  							val |= uint64(input[uint(tableoffset)+bitval*2+1]) << 32
   228  							output[outindex] = val
   229  
   230  							bitpos += uint(encoded_bits)
   231  							outindex++
   232  						}
   233  
   234  					}
   235  				}
   236  			}
   237  
   238  		}
   239  
   240  	}
   241  	return output, nil
   242  }