TWKB编码

"Tiny Well-known Binary" or "TWKB"

Version	Release
0.23	May 1, 2015

Abstract

TWKB is a multi-purpose format for serializing vector geometry data into a byte buffer, with an emphasis on minimizing size of the buffer.

Why not WKB?

The original OGC "well-known binary" format is a simple format, and is capable of easily representing complex OGC geometries like nested collections, but it has two important drawbacks for use as a production serialization:

• it is not aligned, so it doesn't support efficient direct memory access; and,
• it uses IEEE doubles as the coordinate storage format, so for data with lots of spatially adjacent coordinates (basically, all GIS data) it wastes a lot of space on redundant specification of coordinate information.

A new serialization format can address the problem of alignment, or the problem of size, but not both. Given that most current client/server perfomance issues are bottlenecked on network transport times, TWKB concentrates on solving the problem of serialization size.

Basic Principles

TWKB applies the following principles:

• Only store the absolute position once, and store all other positions as delta values relative to the preceding position.
• Only use as much address space as is necessary for any given value. Practically this means that "variable length integers" or "varints" are used throughout the specification for storing values in any situation where numbers greater than 128 might be encountered.

Structure

Standard Attributes

Every TWKB geometry contains standard attributes at the top of the object.

• A type number and precision byte to describe the OGC geometry type.
• A metadata header to indicate which optional attributes to expect, and the storage precision of all coordinates in the geometry.
• An optional extended dimension byte with information about existence and precision of Z & M dimensions.
• An optional size in bytes of the object.
• An optional bounding box of the geometry.
• An optional unique identifier array of sub-components for multi-geometries.

Type & Precision

Size: 1 byte, holding geometry type and number of dimensions

The type-byte stores both the geometry type, and the dimensionality of the coordinates of the geometry.

Bits	Role	Purpose
1-4	Unsigned Integer	Geometry type
5-8	Signed Integer	Geometry precision

• Bits 1-4 store the geometry type (there is space for 16 and we currently use 7):
  • 1 Point
  • 2 Linestring
  • 3 Polygon
  • 4 MultiPoint
  • 5 MultiLinestring
  • 6 MultiPolygon
  • 7 GeometryCollection
• Bits 5-8 store the "precision", which refers to the number of base-10 decimal places stored.
  • A positive precision implies retaining information to the right of the decimal place
    • 41231.1231 at precision=2 would store 41231.12
    • 41231.1231 at precision=1 would store 41231.1
    • 41231.1231 at precision=0 would store 41231
  • A negative precision implies rounding up to the left of the decimal place
    • 41231.1231 at precision=-1 would store 41230
    • 41231.1231 at precision=-2 would store 41200

In order to support negative precisions, the precision number should be stored using zig-zag encoding (see "ZigZag Encode" below).

• The geometry type can be read by masking out the lower four bits (type & 0x0F).
• The precision can be read by masking out the top four bits ((type & 0xF0) >> 4).

Metadata Header

Size: 1 byte

The metadata byte of TWKB is mandatory, and encodes the following information:

Bits	Role	Purpose
1	Boolean	Is there a bounding box?
2	Boolean	Is there a size attribute?
3	Boolean	Is there an ID list?
4	Boolean	Is there extended precision information?
5	Boolean	Is this an empty geometry?
6-8	Boolean	Unused

Extended Dimensions [Optional]

Size: 1 byte

For some coordinate reference systems, and dimension combinations, it makes no sense to store all dimensions using the same precision.

For example, for data with X/Y/Z/T dimensions, using a geographic coordinate system, it might make sense to store the X/Y (longitude/latitude) dimensions with precision of 6 (about 10 cm), the Z with a precision of 1 (also 10 cm), and the T with a precision of 0 (whole seconds). A single precision number cannot handle this case, so the extended precision slot is optionally available for it.

The extended precision byte holds:

Bits	Role	Purpose
1	Boolean	Geometry has Z coordinates?
2	Boolean	Geometry has M coordinates?
3-5	Unsigned Integer	Precision for Z coordinates.
6-8	Unsigned Integer	Precision for M coordinates.

The extended precision values are always positive (only deal with digits to the left of the decimal point)

• The presence of Z can be read by masking out bit 1: (extended & 0x01)
• The presence of M can be read by masking out bit 2: (extended & 0x02)
• The value of Z precision can be read by masking and shifting bits 3-5: (extended & 0x1C) >> 2
• The value of M precision can be read by masking and shifting bits 6-8: (extended & 0xE0) >> 5

Size [Optional]

Size: 1 unsigned varint (so, variable size)

If the size attribute bit is set in the metadata header, a varInt with size information comes next. The value is the size in bytes of the remainder of the geometry after the size attribute.

When encountered in collections, an application can use the size attribute to advance the read pointer the the start of the next geometry. This can be used for a quick scan through a set of geometries, for reading just bounding boxes, or to distibute the read process to different threads.

Bounding Box [Optional]

Size: 2 signed varints per dimension (also variable size)

Each dimension of a bounding box is represented by a pair of varints:

• the minimum of the dimension
• the relative maximum of the dimension (relative to the minimum)

So, for example:

• [xmin, deltax, ymin, deltay, zmin, deltaz]

ID List [Optional]

Size: N signed varints, one per sub-geometry

The TWKB collection types (multipoint, multilinestring, multipolygon, geometrycollection)

• can be used as single values (a given row contains a single collection object), or
• can be used as aggregate wrappers for values from multiple rows (a single object contains geometries read from multiple rows).

In the latter case, it makes sense to include a unique identifier for each sub-geometry that is being wrapped into the collection. The "idlist" attribute is an array of varint that has one varint for each sub-geometry in the collection.

Description of Types

PointArrays

Every object contains a point array, which is an array of coordinates, stored as varints. The final values in the point arrays are a result of four steps:

• convert doubles to integers
• calculate integer delta values between successive points
• zig-zag encode delta values to move negative values into positive range
• varint encode the final value

The result is a very compact representation of the coordinates.

Convert to Integers

All storage using varints, which are integers. However, before being encoded most coordinate are doubles. How are the double coordinates converted into integers? And how are the integers converted ino the final array of varints?

Each coordinate is multiplied by the geometry precision value (from the metadata header), and then rounded to the nearest integer value (round(doublecoord * 10^precision)). When converting from TWKB back to double precision, the reverse process is applied.

Calculate Delta Values

Rather than storing the absolute position of every vertex, we store the difference between each coordinate and the previous coordinate of that dimension. So the coordinate values of a point array are:

x1, y1

(x2 - x1), (y2 - y1)

(x3 - x2), (y3 - y2)

(x4 - x3), (y4 - y3)

...

(xn - xn-1), (yn - yn-1)

Delta values of integer coordinates are nice and small, and they store very compactly as varints.

Every coordinate value in a geometry except the very first one is a delta value relative to the previous value processed. Examples:

• The second coordinate of a line string should be a delta relative to the first.
• The first coordinate of the second ring of a polygon should be a delta relative to the last coordinate of the first ring.
• The first coordinate of the third linestring in a multilinestring should be a delta relative to the last coordinate of the second linestring.

Basically, implementors should always keep the value of the previous coordinate in hand while processing a geometry to use to difference against the incoming coordinate. By setting the very first "previous coordinate" value to [0,0] it is possible to use the same differencing code throughout, since the first value processed will thus receive a "delta" value of itself.

ZigZag Encode

The varint scheme for indicating a value larger than one byte involves setting the "high bit", which is the same thing that standard integers use for storing negative numbers.

So negative numbers need to first be encoded using a "zig zag" encoding, which keep absolute values small, but do not set the high bit as part of encoding negative numbers.

The effect is to convert small negative numbers into small positive numbers, which is exactly what we want, since delta values will tend to be small negative and positive numbers:

-1 => 1

1 => 2

-2 => 3

2 => 4

A computationally fast formula to generate the zig-zag value is

/* for 32-bit signed integer n */

unsigned int zz = (n << 1) ^ (n >> 31)

/* for 64-bit signed integer n */

unsigned long zz = (n << 1) ^ (n >> 63)

VarInt Encode

Variable length integers are a clever way of only using as many bytes as necessary to store a value. The method is described here:

https://developers.google.com/protocol-buffers/docs/encoding#varints

The varint scheme sets the high bit of a byte as a flag to indicate when more bytes are needed to fully represent a number. Decoding a varint accumulates the information in each flagged byte until an unflagged byte is found and the integer is complete and ready to return.

Primitive Types

Point [Type 1]

Because points only have one coordinate, the coordinates will be the true values (not relative) except in cases where the point is embedded in a collection.

Bounding boxes are permitted on points, but discouraged since they just duplicate the values already available in the point coordinate. Similarly, unless the point is part of a collection (where random access is a possible use case), the size should also be omitted.

type_and_prec byte

metadata_header byte

[extended_dims] byte

[size] uvarint

[bounds] bbox

pointarray varint[]

Linestring [Type 2]

A linestring has, in addition to the standard metadata:

• an npoints unsigned varint giving the number of points in the linestring
• if npoints is zero, the linestring is "empty", and there is no further content

The layout is:

type_and_prec byte

metadata_header byte

[extended_dims] byte

[size] uvarint

[bounds] bbox

npoints uvarint

pointarray varint[]

Polygon [Type 3]

A polygon has, in addition to the standard metadata:

• an nrings unsigned varint giving the number of rings in the polygon
• if nrings is zero, the polygon is "empty", and there is no further content
• for each ring there will be
  • an npoints unsigned varint giving the number of points in the ring
  • a pointarray of varints
  • rings are assumed to be implicitly closed, so the first and last point should not be the same

The layout is:

type_and_prec byte

metadata_header byte

[extended_dims] byte

[size] uvarint

[bounds] bbox

nrings uvarint

npoints[0] uvarint

pointarray[0] varint[]

...

npoints[n] uvarint

pointarray[n] varint[]

MultiPoint [Type 4]

A multipoint has, in addition to the standard metadata:

• an optional "idlist" (if indicated in the metadata header)
• an npoints unsigned varint giving the number of points in the multipoint
• if npoints is zero, the multipoint is "empty", and there is no further content

The layout is:

type_and_prec byte

metadata_header byte

[extended_dims] byte

[size] uvarint

[bounds] bbox

npoints uvarint

[idlist] varint[]

pointarray varint[]

MultiLineString [Type 5]

A multilinestring has, in addition to the standard metadata:

• an optional "idlist" (if indicated in the metadata header)
• an nlinestrings unsigned varint giving the number of linestrings in the collection
• if nlinestrings is zero, the collection is "empty", and there is no further content
• for each linestring there will be
  • an npoints unsigned varint giving the number of points in the linestring
  • a pointarray of varints

The layout is:

type_and_prec byte

metadata_header byte

[extended_dims] byte

[size] uvarint

[bounds] bbox

nlinestrings uvarint

[idlist] varint[]

npoints[0] uvarint

pointarray[0] varint[]

...

npoints[n] uvarint

pointarray[n] varint[]

MultiPolygon [Type 6]

A multipolygon has, in addition to the standard metadata:

• an optional "idlist" (if indicated in the metadata header)
• an npolygons unsigned varint giving the number of polygons in the collection
• if npolygons is zero, the collection is "empty", and there is no further content
• for each polygon there will be
  • an nrings unsigned varint giving the number of rings in the linestring
  • for each ring there will be
    • an npoints unsigned varint giving the number of points in the ring
    • a pointarray of varints
    • rings are assumed to be implicitly closed, so the first and last point should not be the same

The layout is:

type_and_prec byte

metadata_header byte

[extended_dims] byte

[size] uvarint

[bounds] bbox

npolygons uvarint

[idlist] varint[]

nrings[0] uvarint

npoints[0][0] uvarint

pointarray[0][0] varint[]

...

nrings[n] uvarint

npoints[n][m] uvarint

pointarray[n][m] varint[]

GeometryCollection [Type 7]

A geometrycollection has, in addition to the standard metadata:

• an optional "idlist" (if indicated in the metadata header)
• an ngeometries unsigned varint giving the number of geometries in the collection
• for each geometry there will be a complete TWKB geometry (with it's own first absolute coordinate), readable using the rules set out above

The layout is:

type_and_prec byte

metadata_header byte

[extended_dims] byte

[size] varint

[bounds] bbox

ngeometries varint

[idlist] varint[]

geom twkb[]

package org.locationtech.geomesa.features.serialization

import com.typesafe.scalalogging.LazyLogging
import org.locationtech.geomesa.utils.geometry.GeometryPrecision.TwkbPrecision
import org.locationtech.jts.geom._

import scala.util.control.NonFatal

/**
  * Based on the TWKB standard: https://github.com/TWKB/Specification/blob/master/twkb.md
  *
  * For backwards compatibility, also reads original serialization, with the `LegacyGeometrySerialization` trait
  */
// noinspection LanguageFeature
trait TwkbSerialization[T <: NumericWriter, V <: NumericReader]
    extends VarIntEncoding[T, V] with WkbSerialization[T, V] with LazyLogging {

  import DimensionalBounds._
  import TwkbSerialization.FlagBytes._
  import TwkbSerialization.GeometryBytes._
  import TwkbSerialization.ZeroByte

  private val factory = new GeometryFactory()
  private val csFactory = factory.getCoordinateSequenceFactory

  /**
    * Serialize a geometry
    *
    * For explanation of precisions, see `org.locationtech.geomesa.utils.geometry.GeometryPrecision`
    *
    * @param out output
    * @param geometry geometry
    * @param precision precision for encoding x, y, z, m
    */
  def serialize(out: T, geometry: Geometry, precision: TwkbPrecision = TwkbPrecision()): Unit = {
    if (geometry == null) {
      out.writeByte(ZeroByte)
    } else {
      // choose our state to correspond with the dimensions in the geometry
      implicit val state: DeltaState = {
        // note that we only check the first coordinate - if a geometry is written with different
        // dimensions in each coordinate, some information may be lost
        val coord = geometry.getCoordinate
        // check for dimensions - use NaN != NaN to verify z coordinate
        // TODO check for M coordinate when added to JTS
        if (coord == null || java.lang.Double.isNaN(coord.getZ)) {
          new XYState(precision.xy)
        } else {
          new XYZState(precision.xy, precision.z)
        }
      }

      geometry match {
        case g: Point =>
          if (g.isEmpty) {
            state.writeMetadata(out, TwkbPoint, empty = true, bbox = false)
          } else {
            state.writeMetadata(out, TwkbPoint, empty = false, bbox = false)
            state.writeCoordinate(out, g.getCoordinate)
          }

        case g: LineString =>
          if (g.isEmpty) {
            state.writeMetadata(out, TwkbLineString, empty = true, bbox = false)
          } else {
            state.writeMetadata(out, TwkbLineString, empty = false, bbox = true)
            state.writeBoundingBox(out, g)
          }
          writeLineString(out, g)

        case g: Polygon =>
          if (g.isEmpty) {
            state.writeMetadata(out, TwkbPolygon, empty = true, bbox = false)
          } else {
            state.writeMetadata(out, TwkbPolygon, empty = false, bbox = true)
            state.writeBoundingBox(out, g)
          }
          writePolygon(out, g)

        case g: MultiPoint =>
          if (g.isEmpty) {
            state.writeMetadata(out, TwkbMultiPoint, empty = true, bbox = false)
          } else {
            state.writeMetadata(out, TwkbMultiPoint, empty = false, bbox = true)
            state.writeBoundingBox(out, g)
          }
          writeMultiPoint(out, g)

        case g: MultiLineString =>
          if (g.isEmpty) {
            state.writeMetadata(out, TwkbMultiLineString, empty = true, bbox = false)
          } else {

            state.writeMetadata(out, TwkbMultiLineString, empty = false, bbox = true)
            state.writeBoundingBox(out, g)
          }
          writeMultiLineString(out, g)

        case g: MultiPolygon =>
          if (g.isEmpty) {
            state.writeMetadata(out, TwkbMultiPolygon, empty = true, bbox = false)
          } else {
            state.writeMetadata(out, TwkbMultiPolygon, empty = false, bbox = true)
            state.writeBoundingBox(out, g)
          }
          writeMultiPolygon(out, g)

        case g: GeometryCollection =>
          if (g.isEmpty) {
            state.writeMetadata(out, TwkbCollection, empty = true, bbox = false)
          } else {
            state.writeMetadata(out, TwkbCollection, empty = false, bbox = true)
            state.writeBoundingBox(out, g)
          }
          writeCollection(out, g)
      }
    }
  }

  /**
    * Deserialize a geometry
    *
    * @param in input
    * @return
    */
  def deserialize(in: V): Geometry = {
    try {
      val precisionAndType = in.readByte()
      if (precisionAndType == ZeroByte) {
        null
      } else if (precisionAndType == NOT_NULL_BYTE) {
        // TODO this overlaps with twkb point type with precision 0
        deserializeWkb(in)
      } else {
        // first byte contains the geometry type in the first 4 bits and the x-y precision in the second 4 bits
        val geomType = (precisionAndType & 0x0F).toByte
        val precision = VarIntEncoding.zigzagDecode((precisionAndType & 0xF0) >>> 4)

        // second byte contains flags for optional elements
        val flags = in.readByte()
        val hasBoundingBox = (flags & BoundingBoxFlag) != 0
        val hasExtendedDims = (flags & ExtendedDimsFlag) != 0
        val isEmpty = (flags & EmptyFlag) != 0

        // extended dims indicates the presence of z and/or m
        // we create our state tracker based on the dimensions that are present
        implicit val state: DeltaState = if (hasExtendedDims) {
          // z and m precisions are indicated in the next byte, where (from right to left):
          //   bit 0 indicates presence of z dimension
          //   bit 1 indicates presence of m dimension
          //   bits 2-5 indicate z precision
          //   bits 6-8 indicate m precision
          val extendedDims = in.readByte()
          if ((extendedDims & 0x01) != 0) { // indicates z dimension
            if ((extendedDims & 0x02) != 0) { // indicates m dimension
              new XYZMState(precision, (extendedDims & 0x1C) >> 2, (extendedDims & 0xE0) >>> 5)
            } else {
              new XYZState(precision, (extendedDims & 0x1C) >> 2)
            }
          } else if ((extendedDims & 0x02) != 0) {  // indicates m dimension
            new XYMState(precision, (extendedDims & 0xE0) >>> 5)
          } else {
            // not sure why anyone would indicate extended dims but set them all false...
            new XYState(precision)
          }
        } else {
          new XYState(precision)
        }

        // size is the length of the remainder of the geometry, after the size attribute
        // we don't currently use size - parsing will fail if size is actually present

        // val hasSize = (flags & FlagBytes.SizeFlag) != 0
        // if (hasSize) {
        //   val size = readUnsignedVarInt(in)
        // }

        // bounding box is not currently used, but we write it in anticipation of future filter optimizations
        if (hasBoundingBox) {
          state.skipBoundingBox(in)
        }

        // children geometries can be written with an id list
        // we don't currently use ids - parsing will fail if ids are actually present
        // val hasIds = (flags & FlagBytes.IdsFlag) != 0

        geomType match {
          case TwkbPoint => factory.createPoint(if (isEmpty) { null } else { csFactory.create(readPointArray(in, 1)) })
          case TwkbLineString      => readLineString(in)
          case TwkbPolygon         => readPolygon(in)
          case TwkbMultiPoint      => readMultiPoint(in)
          case TwkbMultiLineString => readMultiLineString(in)
          case TwkbMultiPolygon    => readMultiPolygon(in)
          case TwkbCollection      => readCollection(in)
          case _ => throw new IllegalArgumentException(s"Invalid TWKB geometry type $geomType")
        }
      }
    } catch {
      case NonFatal(e) => logger.error(s"Error reading serialized kryo geometry:", e); null
    }
  }

  private def writeLineString(out: T, g: LineString)(implicit state: DeltaState): Unit =
      writePointArray(out, g.getCoordinateSequence, g.getNumPoints)

  private def readLineString(in: V)(implicit state: DeltaState): LineString =
    factory.createLineString(csFactory.create(readPointArray(in, readUnsignedVarInt(in))))

  private def writePolygon(out: T, g: Polygon)(implicit state: DeltaState): Unit = {
    if (g.isEmpty) {
      writeUnsignedVarInt(out, 0)
    } else {
      val numRings = g.getNumInteriorRing
      writeUnsignedVarInt(out, numRings + 1) // include exterior ring in count
      // note: don't write final point for each ring, as they should duplicate the first point
      var ring = g.getExteriorRing.getCoordinateSequence
      writePointArray(out, ring, ring.size() - 1)
      var j = 0
      while (j < numRings) {
        ring = g.getInteriorRingN(j).getCoordinateSequence
        writePointArray(out, ring, ring.size() - 1)
        j += 1
      }
    }
  }

  private def readPolygon(in: V)(implicit state: DeltaState): Polygon = {
    val numRings = readUnsignedVarInt(in)
    if (numRings == 0) { factory.createPolygon(null, null) } else {
      val exteriorRing = readLinearRing(in, readUnsignedVarInt(in))
      val interiorRings = Array.ofDim[LinearRing](numRings - 1)
      var i = 1
      while (i < numRings) {
        interiorRings(i - 1) = readLinearRing(in, readUnsignedVarInt(in))
        i += 1
      }
      factory.createPolygon(exteriorRing, interiorRings)
    }
  }

  private def writeMultiPoint(out: T, g: MultiPoint)(implicit state: DeltaState): Unit = {
    val length = g.getNumPoints
    writeUnsignedVarInt(out, length)
    var i = 0
    while (i < length) {
      state.writeCoordinate(out, g.getGeometryN(i).asInstanceOf[Point].getCoordinate)
      i += 1
    }
  }

  private def readMultiPoint(in: V)(implicit state: DeltaState): MultiPoint = {
    val numPoints = readUnsignedVarInt(in)
    if (numPoints == 0) { factory.createMultiPoint(null: CoordinateSequence) } else {
      // note: id list would go here, with one ID per point
      factory.createMultiPoint(readPointArray(in, numPoints).map(factory.createPoint))
    }
  }

  private def writeMultiLineString(out: T, g: MultiLineString)(implicit state: DeltaState): Unit = {
    val length = g.getNumGeometries
    writeUnsignedVarInt(out, length)
    var i = 0
    while (i < length) {
      val line = g.getGeometryN(i).asInstanceOf[LineString].getCoordinateSequence
      writePointArray(out, line, line.size())
      i += 1
    }
  }

  private def readMultiLineString(in: V)(implicit state: DeltaState): MultiLineString = {
    val numLineStrings = readUnsignedVarInt(in)
    if (numLineStrings == 0) { factory.createMultiLineString(null) } else {
      // note: id list would go here, with one ID per line string
      val lineStrings = Array.ofDim[LineString](numLineStrings)
      var i = 0
      while (i < numLineStrings) {
        lineStrings(i) = readLineString(in)
        i += 1
      }
      factory.createMultiLineString(lineStrings)
    }
  }

  private def writeMultiPolygon(out: T, g: MultiPolygon)(implicit state: DeltaState): Unit = {
    val length = g.getNumGeometries
    writeUnsignedVarInt(out, length)
    var i = 0
    while (i < length) {
      writePolygon(out, g.getGeometryN(i).asInstanceOf[Polygon])
      i += 1
    }
  }

  private def readMultiPolygon(in: V)(implicit state: DeltaState): MultiPolygon = {
    val numPolygons = readUnsignedVarInt(in)
    if (numPolygons == 0) { factory.createMultiPolygon(null) } else {
      // note: id list would go here, with one ID per polygon
      val polygons = Array.ofDim[Polygon](numPolygons)
      var i = 0
      while (i < numPolygons) {
        polygons(i) = readPolygon(in)
        i += 1
      }
      factory.createMultiPolygon(polygons)
    }
  }

  private def writeCollection(out: T, g: GeometryCollection)(implicit state: DeltaState): Unit = {
    val length = g.getNumGeometries
    writeUnsignedVarInt(out, length)
    var i = 0
    while (i < length) {
      serialize(out, g.getGeometryN(i))
      i += 1
    }
  }

  private def readCollection(in: V): GeometryCollection = {
    val numGeoms = readUnsignedVarInt(in)
    if (numGeoms == 0) { factory.createGeometryCollection(null) } else {
      // note: id list would go here, with one ID per sub geometry
      val geoms = Array.ofDim[Geometry](numGeoms)
      var i = 0
      while (i < numGeoms) {
        geoms(i) = deserialize(in)
        i += 1
      }
      factory.createGeometryCollection(geoms)
    }
  }

  private def writePointArray(out: T, coords: CoordinateSequence, length: Int)(implicit state: DeltaState): Unit = {
    writeUnsignedVarInt(out, length)
    var i = 0
    while (i < length) {
      state.writeCoordinate(out, coords.getCoordinate(i))
      i += 1
    }
  }

  private def readPointArray(in: V, length: Int)(implicit state: DeltaState): Array[Coordinate] = {
    val result = Array.ofDim[Coordinate](length)
    var i = 0
    while (i < length) {
      result(i) = state.readCoordinate(in)
      i += 1
    }
    result
  }

  private def readLinearRing(in: V, length: Int)(implicit state: DeltaState): LinearRing = {
    if (length == 0) { factory.createLinearRing(null: CoordinateSequence) } else {
      val result = Array.ofDim[Coordinate](length + 1)
      var i = 0
      while (i < length) {
        result(i) = state.readCoordinate(in)
        i += 1
      }
      // linear rings should not store the final, duplicate point, but still need it for the geometry
      result(length) = result(0)
      factory.createLinearRing(csFactory.create(result))
    }
  }

  /**
    * TWKB only reads and writes the delta from one coordinate to the next, which generally saves space
    * over absolute values. This trait tracks the values used for delta calculations over a single
    * read or write operation
    */
  private sealed trait DeltaState {

    /**
      * Write metadata, which includes the geometry and precision byte, the flag byte, and optionally
      * an extended precision byte
      *
      * @param out output
      * @param geometryType geometry type
      * @param empty indicate that the geometry is empty
      * @param bbox indicate that a bbox will be written
      */
    def writeMetadata(out: T, geometryType: Byte, empty: Boolean, bbox: Boolean): Unit

    /**
      * Writes out a bounding box. Each dimension stores a min value and a delta to the max value
      *
      * @param out output
      * @param geometry geometry
      * @param bounds bounds operation
      * @tparam G geometry type
      */
    def writeBoundingBox[G <: Geometry](out: T, geometry: G)(implicit bounds: DimensionalBounds[G]): Unit

    /**
      * Skips over a bounding box. We don't currently use the bounding box when reading
      *
      * @param in in
      */
    def skipBoundingBox(in: V): Unit

    /**
      * Write a coordinate
      *
      * @param out output
      * @param coordinate coordinate
      */
    def writeCoordinate(out: T, coordinate: Coordinate): Unit

    /**
      * Read a coordinate
      *
      * @param in input
      * @return
      */
    def readCoordinate(in: V): Coordinate

    /**
      * Reset the state back to its original state, suitable for re-use
      */
    def reset(): Unit
  }

  private class XYState(precision: Int) extends DeltaState {

    private val p: Double = math.pow(10, precision)
    private var x: Int = 0
    private var y: Int = 0

    protected val boundingBoxFlag: Byte = BoundingBoxFlag
    protected val emptyFlag: Byte = EmptyFlag
    protected val dimensionsFlag: Byte = ZeroByte

    override def writeMetadata(out: T, geometryType: Byte, empty: Boolean, bbox: Boolean): Unit = {
      // write the geometry type and the main precision
      out.writeByte(((VarIntEncoding.zigzagEncode(precision) << 4) | geometryType).toByte)
      // write the flag byte
      out.writeByte(if (bbox) { boundingBoxFlag } else if (empty) { emptyFlag } else { dimensionsFlag })
    }

    override def writeCoordinate(out: T, coordinate: Coordinate): Unit = {
      val cx = math.round(coordinate.x * p).toInt
      val cy = math.round(coordinate.y * p).toInt
      writeVarInt(out, cx - x)
      writeVarInt(out, cy - y)
      x = cx
      y = cy
    }

    override def readCoordinate(in: V): Coordinate = {
      x = x + readVarInt(in)
      y = y + readVarInt(in)
      new Coordinate(x / p, y / p)
    }

    override def writeBoundingBox[G <: Geometry](out: T, geometry: G)(implicit bounds: DimensionalBounds[G]): Unit = {
      val (minX, maxX) = bounds.x(geometry)
      val intX = math.round(minX * p).toInt
      writeVarInt(out, intX)
      writeVarInt(out, math.round(maxX * p).toInt - intX)
      val (minY, maxY) = bounds.y(geometry)
      val intY = math.round(minY * p).toInt
      writeVarInt(out, intY)
      writeVarInt(out, math.round(maxY * p).toInt - intY)
    }

    override def skipBoundingBox(in: V): Unit = {
      skipVarInt(in)
      skipVarInt(in)
      skipVarInt(in)
      skipVarInt(in)
    }

    override def reset(): Unit = {
      x = 0
      y = 0
    }
  }

  private abstract class ExtendedState(precision: Int) extends XYState(precision) {

    protected def extendedDims: Byte

    override protected val boundingBoxFlag: Byte = (ExtendedDimsFlag | BoundingBoxFlag).toByte
    override protected val emptyFlag: Byte = (ExtendedDimsFlag | EmptyFlag).toByte
    override protected val dimensionsFlag: Byte = ExtendedDimsFlag

    override def writeMetadata(out: T, geometryType: Byte, empty: Boolean, bbox: Boolean): Unit = {
      super.writeMetadata(out, geometryType, empty, bbox)
      // write the extended precision values
      out.writeByte(extendedDims)
    }
  }

  private class XYZState(precision: Int, zPrecision: Int) extends ExtendedState(precision) {

    private val pz: Double = math.pow(10, zPrecision)
    private var z: Int = 0

    // sets bits for z dim, and its precisions
    override protected val extendedDims: Byte = (0x01 | ((zPrecision & 0x03) << 2)).toByte

    override def writeBoundingBox[G <: Geometry](out: T, geometry: G)(implicit bounds: DimensionalBounds[G]): Unit = {
      super.writeBoundingBox(out, geometry)
      val (minZ, maxZ) = bounds.z(geometry)
      val intZ = math.round(minZ * pz).toInt
      writeVarInt(out, intZ)
      writeVarInt(out, math.round(maxZ * pz).toInt - intZ)
    }

    override def writeCoordinate(out: T, coordinate: Coordinate): Unit = {
      super.writeCoordinate(out, coordinate)
      val cz = math.round(coordinate.getZ * pz).toInt
      writeVarInt(out, cz - z)
      z = cz
    }

    override def readCoordinate(in: V): Coordinate = {
      val coord = super.readCoordinate(in)
      z = z + readVarInt(in)
      coord.setZ(z / pz)
      coord
    }

    override def skipBoundingBox(in: V): Unit = {
      super.skipBoundingBox(in)
      skipVarInt(in)
      skipVarInt(in)
    }

    override def reset(): Unit = {
      super.reset()
      z = 0
    }
  }

  private class XYMState(precision: Int, mPrecision: Int) extends ExtendedState(precision) {

    private val pm: Double = math.pow(10, mPrecision)
    private var m: Int = 0

    // sets bit for m dim, and its precisions
    override protected val extendedDims: Byte = (0x02 | ((mPrecision & 0x03) << 5)).toByte

    override def writeBoundingBox[G <: Geometry](out: T, geometry: G)(implicit bounds: DimensionalBounds[G]): Unit = {
      super.writeBoundingBox(out, geometry)
      val (minM, maxM) = bounds.m(geometry)
      val intM = math.round(minM * pm).toInt
      writeVarInt(out, intM)
      writeVarInt(out, math.round(maxM * pm).toInt - intM)
    }

    override def writeCoordinate(out: T, coordinate: Coordinate): Unit = {
      super.writeCoordinate(out, coordinate)
      val cm = 0 // TODO math.round(coordinate.m * pm).toInt
      writeVarInt(out, cm - m)
      m = cm
    }

    override def readCoordinate(in: V): Coordinate = {
      val coord = super.readCoordinate(in)
      m = m + readVarInt(in)
      // TODO set m as 4th ordinate when supported by jts
      coord
    }

    override def skipBoundingBox(in: V): Unit = {
      super.skipBoundingBox(in)
      skipVarInt(in)
      skipVarInt(in)
    }

    override def reset(): Unit = {
      super.reset()
      m = 0
    }
  }

  private class XYZMState(precision: Int, zPrecision: Int, mPrecision: Int) extends XYZState(precision, zPrecision) {

    private val pm: Double = math.pow(10, mPrecision)
    private var m: Int = 0

    // sets bits for both z and m dims, and their precisions
    override protected val extendedDims: Byte =
      (0x03 | ((zPrecision & 0x03) << 2) | ((mPrecision & 0x03) << 5)).toByte

    override def writeBoundingBox[G <: Geometry](out: T, geometry: G)(implicit bounds: DimensionalBounds[G]): Unit = {
      super.writeBoundingBox(out, geometry)
      val (minM, maxM) = bounds.m(geometry)
      val intM = math.round(minM * pm).toInt
      writeVarInt(out, intM)
      writeVarInt(out, math.round(maxM * pm).toInt - intM)
    }

    override def writeCoordinate(out: T, coordinate: Coordinate): Unit = {
      super.writeCoordinate(out, coordinate)
      val cm = 0 // TODO math.round(coordinate.m * pm).toInt
      writeVarInt(out, cm - m)
      m = cm
    }

    override def readCoordinate(in: V): Coordinate = {
      val coord = super.readCoordinate(in)
      m = m + readVarInt(in)
      // TODO set m as 4th ordinate when supported by jts
      coord
    }

    override def skipBoundingBox(in: V): Unit = {
      super.skipBoundingBox(in)
      skipVarInt(in)
      skipVarInt(in)
    }

    override def reset(): Unit = {
      super.reset()
      m = 0
    }
  }
}

object TwkbSerialization {

  val MaxPrecision: Byte = 7

  private val ZeroByte: Byte = 0

  // twkb constants
  object GeometryBytes {
    val TwkbPoint           :Byte = 1
    val TwkbLineString      :Byte = 2
    val TwkbPolygon         :Byte = 3
    val TwkbMultiPoint      :Byte = 4
    val TwkbMultiLineString :Byte = 5
    val TwkbMultiPolygon    :Byte = 6
    val TwkbCollection      :Byte = 7
  }

  object FlagBytes {
    val BoundingBoxFlag  :Byte = 0x01
    val SizeFlag         :Byte = 0x02
    val IdsFlag          :Byte = 0x04
    val ExtendedDimsFlag :Byte = 0x08
    val EmptyFlag        :Byte = 0x10
  }
}

-------------------------------------------- C ------------------------------------------

#include "lwout_twkb.h"

/*
* GeometryType, and dimensions
*/
static uint8_t lwgeom_twkb_type(const LWGEOM *geom)
{
    uint8_t twkb_type = 0;

    LWDEBUGF(2, "Entered  lwgeom_twkb_type",0);

    switch ( geom->type )
    {
        case POINTTYPE:
            twkb_type = WKB_POINT_TYPE;
            break;
        case LINETYPE:
            twkb_type = WKB_LINESTRING_TYPE;
            break;
        case TRIANGLETYPE:
        case POLYGONTYPE:
            twkb_type = WKB_POLYGON_TYPE;
            break;
        case MULTIPOINTTYPE:
            twkb_type = WKB_MULTIPOINT_TYPE;
            break;
        case MULTILINETYPE:
            twkb_type = WKB_MULTILINESTRING_TYPE;
            break;
        case MULTIPOLYGONTYPE:
            twkb_type = WKB_MULTIPOLYGON_TYPE;
            break;
        case TINTYPE:
        case COLLECTIONTYPE:
            twkb_type = WKB_GEOMETRYCOLLECTION_TYPE;
            break;
        default:
            lwerror("%s: Unsupported geometry type: %s", __func__, lwtype_name(geom->type));
    }
    return twkb_type;
}


/**
* Calculates the size of the bbox in varints in the form:
* xmin, xdelta, ymin, ydelta
*/
static size_t sizeof_bbox(TWKB_STATE *ts, int ndims)
{
    int i;
    uint8_t buf[16];
    size_t size = 0;
    LWDEBUGF(2, "Entered %s", __func__);
    for ( i = 0; i < ndims; i++ )
    {
        size += varint_s64_encode_buf(ts->bbox_min[i], buf);
        size += varint_s64_encode_buf((ts->bbox_max[i] - ts->bbox_min[i]), buf);
    }
    return size;
}
/**
* Writes the bbox in varints in the form:
* xmin, xdelta, ymin, ydelta
*/
static void write_bbox(TWKB_STATE *ts, int ndims)
{
    int i;
    LWDEBUGF(2, "Entered %s", __func__);
    for ( i = 0; i < ndims; i++ )
    {
        bytebuffer_append_varint(ts->header_buf, ts->bbox_min[i]);
        bytebuffer_append_varint(ts->header_buf, (ts->bbox_max[i] - ts->bbox_min[i]));
    }
}


/**
* Stores a pointarray as varints in the buffer
* @register_npoints, controls whether an npoints entry is added to the buffer (used to skip npoints for point types)
* @dimension, states the dimensionality of object this array is part of (0 = point, 1 = linear, 2 = areal)
*/
static int ptarray_to_twkb_buf(const POINTARRAY *pa, TWKB_GLOBALS *globals, TWKB_STATE *ts, int register_npoints, uint32_t minpoints)
{
    uint32_t ndims = FLAGS_NDIMS(pa->flags);
    uint32_t i, j;
    bytebuffer_t b;
    bytebuffer_t *b_p;
    int64_t nextdelta[MAX_N_DIMS];
    int npoints = 0;
    size_t npoints_offset = 0;
    uint32_t max_points_left = pa->npoints;

    LWDEBUGF(2, "Entered %s", __func__);

    /* Dispense with the empty case right away */
    if ( pa->npoints == 0 && register_npoints )
    {
        LWDEBUGF(4, "Register npoints:%d", pa->npoints);
        bytebuffer_append_uvarint(ts->geom_buf, pa->npoints);
        return 0;
    }

    /* If npoints is more than 127 it is unpredictable how many bytes npoints will need */
    /* Then we have to store the deltas in a temp buffer to later add them after npoints */
    /* If noints is below 128 we know 1 byte will be needed */
    /* Then we can make room for that 1 byte at once and write to */
    /* ordinary buffer */
    if( pa->npoints > 127 )
    {
        /* Independent buffer to hold the coordinates, so we can put the npoints */
        /* into the stream once we know how many points we actually have */
        bytebuffer_init_with_size(&b, 3 * ndims * pa->npoints);
        b_p = &b;
    }
    else
    {
        /* We give an alias to our ordinary buffer */
        b_p = ts->geom_buf;
        if ( register_npoints )
        {
            /* We do not store a pointer to the place where we want the npoints value */
            /* Instead we store how far from the beginning of the buffer we want the value */
            /* That is because we otherwise will get in trouble if the buffer is reallocated */
            npoints_offset = b_p->writecursor - b_p->buf_start;

            /* We just move the cursor 1 step to make room for npoints byte */
            /* We use the function append_byte even if we have no value yet, */
            /* since that gives us the check for big enough buffer and moves the cursor */
            bytebuffer_append_byte(b_p, 0);
        }
    }

    for ( i = 0; i < pa->npoints; i++ )
    {
        double *dbl_ptr = (double*)getPoint_internal(pa, i);
        int64_t diff = 0;

        /* Write this coordinate to the buffer as a varint */
        for ( j = 0; j < ndims; j++ )
        {
            /* To get the relative coordinate we don't get the distance */
            /* from the last point but instead the distance from our */
            /* last accumulated point. This is important to not build up an */
            /* accumulated error when rounding the coordinates */
            nextdelta[j] = (int64_t) llround(globals->factor[j] * dbl_ptr[j]) - ts->accum_rels[j];
            LWDEBUGF(4, "deltavalue: %d, ", nextdelta[j]);
            diff += llabs(nextdelta[j]);
        }

        /* Skipping the first point is not allowed */
        /* If the sum(abs()) of all the deltas was zero, */
        /* then this was a duplicate point, so we can ignore it */
        if ( i > 0 && diff == 0 &&  max_points_left > minpoints )
        {
            max_points_left--;
            continue;
        }

        /* We really added a point, so... */
        npoints++;

        /* Write this vertex to the temporary buffer as varints */
        for ( j = 0; j < ndims; j++ )
        {
            ts->accum_rels[j] += nextdelta[j];
            bytebuffer_append_varint(b_p, nextdelta[j]);
        }

        /* See if this coordinate expands the bounding box */
        if( globals->variant & TWKB_BBOX )
        {
            for ( j = 0; j < ndims; j++ )
            {
                if( ts->accum_rels[j] > ts->bbox_max[j] )
                    ts->bbox_max[j] = ts->accum_rels[j];

                if( ts->accum_rels[j] < ts->bbox_min[j] )
                    ts->bbox_min[j] = ts->accum_rels[j];
            }
        }

    }

    if ( pa->npoints > 127 )
    {
        /* Now write the temporary results into the main buffer */
        /* First the npoints */
        if ( register_npoints )
            bytebuffer_append_uvarint(ts->geom_buf, npoints);
        /* Now the coordinates */
        bytebuffer_append_bytebuffer(ts->geom_buf, b_p);

        /* Clear our temporary buffer */
        bytebuffer_destroy_buffer(&b);
    }
    else
    {
        /* If we didn't use a temp buffer, we just write that npoints value */
        /* to where it belongs*/
        if ( register_npoints )
            varint_u64_encode_buf(npoints, b_p->buf_start + npoints_offset);
    }

    return 0;
}

/******************************************************************
* POINTS
*******************************************************************/

static int lwpoint_to_twkb_buf(const LWPOINT *pt, TWKB_GLOBALS *globals, TWKB_STATE *ts)
{
    LWDEBUGF(2, "Entered %s", __func__);

    /* Set the coordinates (don't write npoints) */
    ptarray_to_twkb_buf(pt->point, globals, ts, 0, 1);
    return 0;
}

/******************************************************************
* LINESTRINGS
*******************************************************************/

static int lwline_to_twkb_buf(const LWLINE *line, TWKB_GLOBALS *globals, TWKB_STATE *ts)
{
    LWDEBUGF(2, "Entered %s", __func__);

    /* Set the coordinates (do write npoints) */
    ptarray_to_twkb_buf(line->points, globals, ts, 1, 2);
    return 0;
}

static int
lwtriangle_to_twkb_buf(const LWTRIANGLE *tri, TWKB_GLOBALS *globals, TWKB_STATE *ts)
{
    LWDEBUGF(2, "Entered %s", __func__);
    bytebuffer_append_uvarint(ts->geom_buf, (uint64_t)1);

    /* Set the coordinates (do write npoints) */
    ptarray_to_twkb_buf(tri->points, globals, ts, 1, 2);
    return 0;
}

/******************************************************************
* POLYGONS
*******************************************************************/

static int lwpoly_to_twkb_buf(const LWPOLY *poly, TWKB_GLOBALS *globals, TWKB_STATE *ts)
{
    uint32_t i;

    /* Set the number of rings */
    bytebuffer_append_uvarint(ts->geom_buf, (uint64_t) poly->nrings);

    for ( i = 0; i < poly->nrings; i++ )
    {
        /* Set the coordinates (do write npoints) */
        ptarray_to_twkb_buf(poly->rings[i], globals, ts, 1, 4);
    }

    return 0;
}



/******************************************************************
* MULTI-GEOMETRYS (MultiPoint, MultiLinestring, MultiPolygon)
*******************************************************************/

static int lwmulti_to_twkb_buf(const LWCOLLECTION *col, TWKB_GLOBALS *globals, TWKB_STATE *ts)
{
    uint32_t i;
    int nempty = 0;

    LWDEBUGF(2, "Entered %s", __func__);
    LWDEBUGF(4, "Number of geometries in multi is %d", col->ngeoms);

    /* Deal with special case for MULTIPOINT: skip any empty points */
    if ( col->type == MULTIPOINTTYPE )
    {
        for ( i = 0; i < col->ngeoms; i++ )
            if ( lwgeom_is_empty(col->geoms[i]) )
                nempty++;
    }

    /* Set the number of geometries */
    bytebuffer_append_uvarint(ts->geom_buf, (uint64_t) (col->ngeoms - nempty));

    /* We've been handed an idlist, so write it in */
    if ( ts->idlist )
    {
        for ( i = 0; i < col->ngeoms; i++ )
        {
            /* Skip empty points in multipoints, we can't represent them */
            if ( col->type == MULTIPOINTTYPE && lwgeom_is_empty(col->geoms[i]) )
                continue;

            bytebuffer_append_varint(ts->geom_buf, ts->idlist[i]);
        }

        /* Empty it out to nobody else uses it now */
        ts->idlist = NULL;
    }

    for ( i = 0; i < col->ngeoms; i++ )
    {
        /* Skip empty points in multipoints, we can't represent them */
        if ( col->type == MULTIPOINTTYPE && lwgeom_is_empty(col->geoms[i]) )
            continue;

        lwgeom_to_twkb_buf(col->geoms[i], globals, ts);
    }
    return 0;
}

/******************************************************************
* GEOMETRYCOLLECTIONS
*******************************************************************/

static int lwcollection_to_twkb_buf(const LWCOLLECTION *col, TWKB_GLOBALS *globals, TWKB_STATE *ts)
{
    uint32_t i;

    LWDEBUGF(2, "Entered %s", __func__);
    LWDEBUGF(4, "Number of geometries in collection is %d", col->ngeoms);

    /* Set the number of geometries */
    bytebuffer_append_uvarint(ts->geom_buf, (uint64_t) col->ngeoms);

    /* We've been handed an idlist, so write it in */
    if ( ts->idlist )
    {
        for ( i = 0; i < col->ngeoms; i++ )
            bytebuffer_append_varint(ts->geom_buf, ts->idlist[i]);

        /* Empty it out to nobody else uses it now */
        ts->idlist = NULL;
    }

    /* Write in the sub-geometries */
    for ( i = 0; i < col->ngeoms; i++ )
    {
        lwgeom_write_to_buffer(col->geoms[i], globals, ts);
    }
    return 0;
}


/******************************************************************
* Handle whole TWKB
*******************************************************************/

static int lwgeom_to_twkb_buf(const LWGEOM *geom, TWKB_GLOBALS *globals, TWKB_STATE *ts)
{
    LWDEBUGF(2, "Entered %s", __func__);

    switch ( geom->type )
    {
        case POINTTYPE:
        {
            LWDEBUGF(4,"Type found is Point, %d", geom->type);
            return lwpoint_to_twkb_buf((LWPOINT*) geom, globals, ts);
        }
        case LINETYPE:
        {
            LWDEBUGF(4,"Type found is Linestring, %d", geom->type);
            return lwline_to_twkb_buf((LWLINE*) geom, globals, ts);
        }
        case TRIANGLETYPE:
        {
            LWDEBUGF(4, "Type found is Triangle, %d", geom->type);
            return lwtriangle_to_twkb_buf((LWTRIANGLE *)geom, globals, ts);
        }
        /* Polygon has 'nrings' and 'rings' elements */
        case POLYGONTYPE:
        {
            LWDEBUGF(4,"Type found is Polygon, %d", geom->type);
            return lwpoly_to_twkb_buf((LWPOLY*)geom, globals, ts);
        }

        /* All these Collection types have 'ngeoms' and 'geoms' elements */
        case MULTIPOINTTYPE:
        case MULTILINETYPE:
        case MULTIPOLYGONTYPE:
        {
            LWDEBUGF(4,"Type found is Multi, %d", geom->type);
            return lwmulti_to_twkb_buf((LWCOLLECTION*)geom, globals, ts);
        }
        case COLLECTIONTYPE:
        case TINTYPE:
        {
            LWDEBUGF(4,"Type found is collection, %d", geom->type);
            return lwcollection_to_twkb_buf((LWCOLLECTION*) geom, globals, ts);
        }
        /* Unknown type! */
        default:
            lwerror("%s: Unsupported geometry type: %s", __func__, lwtype_name(geom->type));
    }

    return 0;
}


static int lwgeom_write_to_buffer(const LWGEOM *geom, TWKB_GLOBALS *globals, TWKB_STATE *parent_state)
{
    int i, is_empty, has_z = 0, has_m = 0, ndims;
    size_t bbox_size = 0, optional_precision_byte = 0;
    uint8_t flag = 0, type_prec = 0;
    bytebuffer_t header_bytebuffer, geom_bytebuffer;

    TWKB_STATE child_state;
    memset(&child_state, 0, sizeof(TWKB_STATE));
    child_state.header_buf = &header_bytebuffer;
    child_state.geom_buf = &geom_bytebuffer;
    child_state.idlist = parent_state->idlist;

    bytebuffer_init_with_size(child_state.header_buf, 16);
    bytebuffer_init_with_size(child_state.geom_buf, 64);

    /* Read dimensionality from input */
    ndims = lwgeom_ndims(geom);
    is_empty = lwgeom_is_empty(geom);
    if ( ndims > 2 )
    {
        has_z = lwgeom_has_z(geom);
        has_m = lwgeom_has_m(geom);
    }

    /* Do we need extended precision? If we have a Z or M we do. */
    optional_precision_byte = (has_z || has_m);

    /* Both X and Y dimension use the same precision */
    globals->factor[0] = pow(10, globals->prec_xy);
    globals->factor[1] = globals->factor[0];

    /* Z and M dimensions have their own precisions */
    if ( has_z )
        globals->factor[2] = pow(10, globals->prec_z);
    if ( has_m )
        globals->factor[2 + has_z] = pow(10, globals->prec_m);

    /* Reset stats */
    for ( i = 0; i < MAX_N_DIMS; i++ )
    {
        /* Reset bbox calculation */
        child_state.bbox_max[i] = INT64_MIN;
        child_state.bbox_min[i] = INT64_MAX;
        /* Reset acumulated delta values to get absolute values on next point */
        child_state.accum_rels[i] = 0;
    }

    /* TYPE/PRECISION BYTE */
    if ( abs(globals->prec_xy) > 7 )
        lwerror("%s: X/Z precision cannot be greater than 7 or less than -7", __func__);

    /* Read the TWKB type number from the geometry */
    TYPE_PREC_SET_TYPE(type_prec, lwgeom_twkb_type(geom));
    /* Zig-zag the precision value before encoding it since it is a signed value */
    TYPE_PREC_SET_PREC(type_prec, zigzag8(globals->prec_xy));
    /* Write the type and precision byte */
    bytebuffer_append_byte(child_state.header_buf, type_prec);

    /* METADATA BYTE */
    /* Set first bit if we are going to store bboxes */
    FIRST_BYTE_SET_BBOXES(flag, (globals->variant & TWKB_BBOX) && ! is_empty);
    /* Set second bit if we are going to store resulting size */
    FIRST_BYTE_SET_SIZES(flag, globals->variant & TWKB_SIZE);
    /* There will be no ID-list (for now) */
    FIRST_BYTE_SET_IDLIST(flag, parent_state->idlist && ! is_empty);
    /* Are there higher dimensions */
    FIRST_BYTE_SET_EXTENDED(flag, optional_precision_byte);
    /* Empty? */
    FIRST_BYTE_SET_EMPTY(flag, is_empty);
    /* Write the header byte */
    bytebuffer_append_byte(child_state.header_buf, flag);

    /* EXTENDED PRECISION BYTE (OPTIONAL) */
    /* If needed, write the extended dim byte */
    if( optional_precision_byte )
    {
        uint8_t flag = 0;

        if ( has_z && ( globals->prec_z > 7 || globals->prec_z < 0 ) )
            lwerror("%s: Z precision cannot be negative or greater than 7", __func__);

        if ( has_m && ( globals->prec_m > 7 || globals->prec_m < 0 ) )
            lwerror("%s: M precision cannot be negative or greater than 7", __func__);

        HIGHER_DIM_SET_HASZ(flag, has_z);
        HIGHER_DIM_SET_HASM(flag, has_m);
        HIGHER_DIM_SET_PRECZ(flag, globals->prec_z);
        HIGHER_DIM_SET_PRECM(flag, globals->prec_m);
        bytebuffer_append_byte(child_state.header_buf, flag);
    }

    /* It the geometry is empty, we're almost done */
    if ( is_empty )
    {
        /* If this output is sized, write the size of */
        /* all following content, which is zero because */
        /* there is none */
        if ( globals->variant & TWKB_SIZE )
            bytebuffer_append_byte(child_state.header_buf, 0);

        bytebuffer_append_bytebuffer(parent_state->geom_buf, child_state.header_buf);
        bytebuffer_destroy_buffer(child_state.header_buf);
        bytebuffer_destroy_buffer(child_state.geom_buf);
        return 0;
    }

    /* Write the TWKB into the output buffer */
    lwgeom_to_twkb_buf(geom, globals, &child_state);

    /*If we have a header_buf, we know that this function is called inside a collection*/
    /*and then we have to merge the bboxes of the included geometries*/
    /*and put the result to the parent (the collection)*/
    if( (globals->variant & TWKB_BBOX) && parent_state->header_buf )
    {
        LWDEBUG(4,"Merge bboxes");
        for ( i = 0; i < ndims; i++ )
        {
            if(child_state.bbox_min[i]<parent_state->bbox_min[i])
                parent_state->bbox_min[i] = child_state.bbox_min[i];
            if(child_state.bbox_max[i]>parent_state->bbox_max[i])
                parent_state->bbox_max[i] = child_state.bbox_max[i];
        }
    }

    /* Did we have a box? If so, how big? */
    bbox_size = 0;
    if( globals->variant & TWKB_BBOX )
    {
        LWDEBUG(4,"We want boxes and will calculate required size");
        bbox_size = sizeof_bbox(&child_state, ndims);
    }

    /* Write the size if wanted */
    if( globals->variant & TWKB_SIZE )
    {
        /* Here we have to add what we know will be written to header */
        /* buffer after size value is written */
        size_t size_to_register = bytebuffer_getlength(child_state.geom_buf);
        size_to_register += bbox_size;
        bytebuffer_append_uvarint(child_state.header_buf, size_to_register);
    }

    if( globals->variant & TWKB_BBOX )
        write_bbox(&child_state, ndims);

    bytebuffer_append_bytebuffer(parent_state->geom_buf,child_state.header_buf);
    bytebuffer_append_bytebuffer(parent_state->geom_buf,child_state.geom_buf);

    bytebuffer_destroy_buffer(child_state.header_buf);
    bytebuffer_destroy_buffer(child_state.geom_buf);
    return 0;
}


/**
* Convert LWGEOM to a char* in TWKB format. Caller is responsible for freeing
* the returned array.
*/
lwvarlena_t *
lwgeom_to_twkb_with_idlist(const LWGEOM *geom,
               int64_t *idlist,
               uint8_t variant,
               int8_t precision_xy,
               int8_t precision_z,
               int8_t precision_m)
{
    LWDEBUGF(2, "Entered %s", __func__);
    LWDEBUGF(2, "variant value %x", variant);

    TWKB_GLOBALS tg;
    TWKB_STATE ts;
    bytebuffer_t geom_bytebuffer;

    memset(&ts, 0, sizeof(TWKB_STATE));
    memset(&tg, 0, sizeof(TWKB_GLOBALS));

    tg.variant = variant;
    tg.prec_xy = precision_xy;
    tg.prec_z = precision_z;
    tg.prec_m = precision_m;

    if ( idlist && ! lwgeom_is_collection(geom) )
    {
        lwerror("Only collections can support ID lists");
        return NULL;
    }

    if ( ! geom )
    {
        LWDEBUG(4,"Cannot convert NULL into TWKB.");
        lwerror("Cannot convert NULL into TWKB");
        return NULL;
    }

    ts.idlist = idlist;
    ts.header_buf = NULL;
    ts.geom_buf = &geom_bytebuffer;
    bytebuffer_init_with_size(ts.geom_buf, 512);
    lwgeom_write_to_buffer(geom, &tg, &ts);

    lwvarlena_t *v = bytebuffer_get_buffer_varlena(ts.geom_buf);
    bytebuffer_destroy_buffer(ts.geom_buf);
    return v;
}

lwvarlena_t *
lwgeom_to_twkb(const LWGEOM *geom, uint8_t variant, int8_t precision_xy, int8_t precision_z, int8_t precision_m)
{
    return lwgeom_to_twkb_with_idlist(geom, NULL, variant, precision_xy, precision_z, precision_m);
}