From 43d5e17d5a2a18deae60ca1d0c8e34288b0f1d3f Mon Sep 17 00:00:00 2001 From: Thomas Bonfort Date: Tue, 28 Sep 2021 15:57:28 +0200 Subject: [PATCH] support planarconfig=2, fix tile ordering issues (#3) * support planarconfig=2, fix tile ordering issues - support for band interleaved sources / output - fix unoptimized tile interleaving when overviews are present - relax sanity check on geotransform scale --- mucog.go | 242 +++++++++++++++++++++++++---------------------- testdata/main.go | 73 ++++++++++++++ 2 files changed, 202 insertions(+), 113 deletions(-) create mode 100644 testdata/main.go diff --git a/mucog.go b/mucog.go index 93be354..a08716e 100644 --- a/mucog.go +++ b/mucog.go @@ -114,6 +114,7 @@ type IFD struct { ntags uint64 tagsSize uint64 strileSize uint64 + nplanes uint64 //1 if PlanarConfiguration==1, SamplesPerPixel if PlanarConfiguration==2 ntilesx, ntilesy uint64 minx, miny, maxx, maxy uint64 r tiff.BReader @@ -149,84 +150,88 @@ func (ifd *IFD) AddOverview(ovr *IFD) { ifd.SubIFDs = append(ifd.SubIFDs, ovr) } -func (ifd *IFD) structure(bigtiff bool) (tagCount, ifdSize, strileSize uint64) { - cnt := uint64(0) - size := uint64(16) //8 for field count + 8 for next ifd offset +func (ifd *IFD) structure(bigtiff bool) (tagCount, ifdSize, strileSize, planeCount uint64) { + tagCount = 0 + ifdSize = 16 //8 for field count + 8 for next ifd offset tagSize := uint64(20) + planeCount = 1 if !bigtiff { - size = 6 // 2 for field count + 4 for next ifd offset + ifdSize = 6 // 2 for field count + 4 for next ifd offset tagSize = 12 } strileSize = uint64(0) if ifd.SubfileType > 0 { - cnt++ - size += tagSize + tagCount++ + ifdSize += tagSize } if ifd.ImageWidth > 0 { - cnt++ - size += tagSize + tagCount++ + ifdSize += tagSize } if ifd.ImageLength > 0 { - cnt++ - size += tagSize + tagCount++ + ifdSize += tagSize } if len(ifd.BitsPerSample) > 0 { - cnt++ - size += arrayFieldSize(ifd.BitsPerSample, bigtiff) + tagCount++ + ifdSize += arrayFieldSize(ifd.BitsPerSample, bigtiff) } if ifd.Compression > 0 { - cnt++ - size += tagSize + tagCount++ + ifdSize += tagSize } - cnt++ /*PhotometricInterpretation*/ - size += tagSize + tagCount++ /*PhotometricInterpretation*/ + ifdSize += tagSize if len(ifd.DocumentName) > 0 { - cnt++ - size += arrayFieldSize(ifd.DocumentName, bigtiff) + tagCount++ + ifdSize += arrayFieldSize(ifd.DocumentName, bigtiff) } if ifd.SamplesPerPixel > 0 { - cnt++ - size += tagSize + tagCount++ + ifdSize += tagSize } if ifd.PlanarConfiguration > 0 { - cnt++ - size += tagSize + tagCount++ + ifdSize += tagSize + if ifd.PlanarConfiguration == PlanarConfigurationSeparate { + planeCount = uint64(ifd.SamplesPerPixel) + } } if len(ifd.DateTime) > 0 { - cnt++ - size += arrayFieldSize(ifd.DateTime, bigtiff) + tagCount++ + ifdSize += arrayFieldSize(ifd.DateTime, bigtiff) } if ifd.Predictor > 0 { - cnt++ - size += tagSize + tagCount++ + ifdSize += tagSize } if len(ifd.Colormap) > 0 { - cnt++ - size += arrayFieldSize(ifd.BitsPerSample, bigtiff) + tagCount++ + ifdSize += arrayFieldSize(ifd.BitsPerSample, bigtiff) } if ifd.TileWidth > 0 { - cnt++ - size += tagSize + tagCount++ + ifdSize += tagSize } if ifd.TileLength > 0 { - cnt++ - size += tagSize + tagCount++ + ifdSize += tagSize } if len(ifd.NewTileOffsets32) > 0 { - cnt++ - size += tagSize + tagCount++ + ifdSize += tagSize strileSize += arrayFieldSize(ifd.NewTileOffsets32, bigtiff) - tagSize } else if len(ifd.NewTileOffsets64) > 0 { - cnt++ - size += tagSize + tagCount++ + ifdSize += tagSize strileSize += arrayFieldSize(ifd.NewTileOffsets64, bigtiff) - tagSize } if len(ifd.TileByteCounts) > 0 { - cnt++ - size += tagSize + tagCount++ + ifdSize += tagSize strileSize += arrayFieldSize(ifd.TileByteCounts, bigtiff) - tagSize } if len(ifd.SubIFDOffsets) > 0 { @@ -234,62 +239,62 @@ func (ifd *IFD) structure(bigtiff bool) (tagCount, ifdSize, strileSize uint64) { for i := range offs { offs[i] = uint32(ifd.SubIFDOffsets[i]) } - cnt++ - size += arrayFieldSize(offs, bigtiff) + tagCount++ + ifdSize += arrayFieldSize(offs, bigtiff) } if len(ifd.ExtraSamples) > 0 { - cnt++ - size += arrayFieldSize(ifd.ExtraSamples, bigtiff) + tagCount++ + ifdSize += arrayFieldSize(ifd.ExtraSamples, bigtiff) } if len(ifd.SampleFormat) > 0 { - cnt++ - size += arrayFieldSize(ifd.SampleFormat, bigtiff) + tagCount++ + ifdSize += arrayFieldSize(ifd.SampleFormat, bigtiff) } if len(ifd.JPEGTables) > 0 { - cnt++ - size += arrayFieldSize(ifd.JPEGTables, bigtiff) + tagCount++ + ifdSize += arrayFieldSize(ifd.JPEGTables, bigtiff) } if len(ifd.ModelPixelScaleTag) > 0 { - cnt++ - size += arrayFieldSize(ifd.ModelPixelScaleTag, bigtiff) + tagCount++ + ifdSize += arrayFieldSize(ifd.ModelPixelScaleTag, bigtiff) } if len(ifd.ModelTiePointTag) > 0 { - cnt++ - size += arrayFieldSize(ifd.ModelTiePointTag, bigtiff) + tagCount++ + ifdSize += arrayFieldSize(ifd.ModelTiePointTag, bigtiff) } if len(ifd.ModelTransformationTag) > 0 { - cnt++ - size += arrayFieldSize(ifd.ModelTransformationTag, bigtiff) + tagCount++ + ifdSize += arrayFieldSize(ifd.ModelTransformationTag, bigtiff) } if len(ifd.GeoKeyDirectoryTag) > 0 { - cnt++ - size += arrayFieldSize(ifd.GeoKeyDirectoryTag, bigtiff) + tagCount++ + ifdSize += arrayFieldSize(ifd.GeoKeyDirectoryTag, bigtiff) } if len(ifd.GeoDoubleParamsTag) > 0 { - cnt++ - size += arrayFieldSize(ifd.GeoDoubleParamsTag, bigtiff) + tagCount++ + ifdSize += arrayFieldSize(ifd.GeoDoubleParamsTag, bigtiff) } if ifd.GeoAsciiParamsTag != "" { - cnt++ - size += arrayFieldSize(ifd.GeoAsciiParamsTag, bigtiff) + tagCount++ + ifdSize += arrayFieldSize(ifd.GeoAsciiParamsTag, bigtiff) } if ifd.GDALMetaData != "" { - cnt++ - size += arrayFieldSize(ifd.GDALMetaData, bigtiff) + tagCount++ + ifdSize += arrayFieldSize(ifd.GDALMetaData, bigtiff) } if len(ifd.LERCParams) > 0 { - cnt++ - size += arrayFieldSize(ifd.LERCParams, bigtiff) + tagCount++ + ifdSize += arrayFieldSize(ifd.LERCParams, bigtiff) } if len(ifd.RPCs) > 0 { - cnt++ - size += arrayFieldSize(ifd.RPCs, bigtiff) + tagCount++ + ifdSize += arrayFieldSize(ifd.RPCs, bigtiff) } if ifd.NoData != "" { - cnt++ - size += arrayFieldSize(ifd.NoData, bigtiff) + tagCount++ + ifdSize += arrayFieldSize(ifd.NoData, bigtiff) } - return cnt, size, strileSize + return } type TagData struct { @@ -384,18 +389,23 @@ func (cog *MultiCOG) computeStructure(bigtiff bool) error { */ for i, ifd := range cog.ifds { - ifd.ntags, ifd.tagsSize, ifd.strileSize = ifd.structure(bigtiff) + ifd.ntags, ifd.tagsSize, ifd.strileSize, ifd.nplanes = ifd.structure(bigtiff) ifd.ntilesx = (ifd.ImageWidth + uint64(ifd.TileWidth) - 1) / uint64(ifd.TileWidth) ifd.ntilesy = (ifd.ImageLength + uint64(ifd.TileLength) - 1) / uint64(ifd.TileLength) isx, isy := ifd.gt.Scale() - if isx != sx || isy != sy { - return fmt.Errorf("ifd %d incompatible scale (x: %f/%f, y: %f/%f)", i, isx, sx, isy, sy) + xScaleDiff := math.Abs(1 - isx/sx) + yScaleDiff := math.Abs(1 - isy/sy) + if xScaleDiff > 0.00000001 || yScaleDiff > 0.00000001 { + return fmt.Errorf("ifd %d incompatible scales (x: %.16f/%.16f, y: %.16f/%.16f)", i, isx, sx, isy, sy) } if ifd.TileWidth != tsx || ifd.TileLength != tsy { return fmt.Errorf("ifd %d incompatible tile size (sx: %d/%d, sy: %d/%d)", i, ifd.TileWidth, tsx, ifd.TileLength, tsy) } + if ifd.nplanes != cog.ifds[0].nplanes { + return fmt.Errorf("ifd %d incompatible number of planes (%d/%d)", i, ifd.nplanes, cog.ifds[0].nplanes) + } iox, ioy := ifd.gt.Origin() //pixel offset from origin of first ifd @@ -412,11 +422,11 @@ func (cog *MultiCOG) computeStructure(bigtiff bool) error { ifd.maxx = ifd.minx + ifd.ntilesx ifd.maxy = ifd.miny + ifd.ntilesy - for _, ifd := range ifd.SubIFDs { - ifd.ntags, ifd.tagsSize, ifd.strileSize = ifd.structure(bigtiff) - ifd.ntilesx = (ifd.ImageWidth + uint64(ifd.TileWidth) - 1) / uint64(ifd.TileWidth) - ifd.ntilesy = (ifd.ImageLength + uint64(ifd.TileLength) - 1) / uint64(ifd.TileLength) - ifd.minx, ifd.miny, ifd.maxx, ifd.maxy = 0, 0, ifd.ntilesx, ifd.ntilesy + for _, sifd := range ifd.SubIFDs { + sifd.ntags, sifd.tagsSize, sifd.strileSize, sifd.nplanes = sifd.structure(bigtiff) + sifd.ntilesx = (sifd.ImageWidth + uint64(sifd.TileWidth) - 1) / uint64(sifd.TileWidth) + sifd.ntilesy = (sifd.ImageLength + uint64(sifd.TileLength) - 1) / uint64(sifd.TileLength) + sifd.minx, sifd.miny, sifd.maxx, sifd.maxy = 0, 0, sifd.ntilesx, sifd.ntilesy } } return nil @@ -466,7 +476,7 @@ func (cog *MultiCOG) computeImageryOffsets(bigtiff bool) error { datas := cog.dataInterlacing() tiles := datas.Tiles() for tile := range tiles { - tileidx := tile.x + tile.y*tile.ifd.ntilesx + tileidx := (tile.x+tile.y*tile.ifd.ntilesx)*tile.ifd.nplanes + tile.plane cnt := uint64(tile.ifd.TileByteCounts[tileidx]) if cnt > 0 { if bigtiff { @@ -557,7 +567,7 @@ func (cog *MultiCOG) Write(out io.Writer, bigtiff bool) error { buf := &bytes.Buffer{} for tile := range tiles { buf.Reset() - idx := tile.x + tile.y*tile.ifd.ntilesx + idx := (tile.x+tile.y*tile.ifd.ntilesx)*tile.ifd.nplanes + tile.plane if tile.ifd.TileByteCounts[idx] > 0 { _, err := tile.ifd.r.Seek(int64(tile.ifd.OriginalTileOffsets[idx]), io.SeekStart) if err != nil { @@ -847,19 +857,22 @@ func (cog *MultiCOG) writeIFD(w io.Writer, bigtiff bool, ifd *IFD, offset uint64 } type tile struct { - ifd *IFD - x, y uint64 + ifd *IFD + x, y uint64 + plane uint64 } func (cog *MultiCOG) dataInterlacing() datas { - ret := make([][]*IFD, 1) + ret := [][]*IFD{} + full := []*IFD{} + for _, topifd := range cog.ifds { - ret[0] = append(ret[0], topifd) + full = append(full, topifd) ovr := []*IFD{} for _, subifd := range topifd.SubIFDs { if subifd.SubfileType == SubfileTypeMask && subifd.ImageWidth == topifd.ImageWidth { - ret[0] = append(ret[0], subifd) + full = append(full, subifd) } else { ovr = append(ovr, subifd) } @@ -874,9 +887,7 @@ func (cog *MultiCOG) dataInterlacing() datas { ret = append(ret, ovr) } } - //send main ifd to the end - ret = append(ret, ret[0]) - ret = ret[1:] + ret = append(ret, full) return ret } @@ -887,50 +898,55 @@ func (d datas) Tiles() chan tile { go func() { defer close(ch) - ifds := d[0] - y := uint64(0) - for { - yok := false - for _, ifd := range ifds { - if y < ifd.maxy { - yok = true - break - } - } - if !yok { - break - } - x := uint64(0) + ifds := d[len(d)-1] + for plane := uint64(0); plane < ifds[0].nplanes; plane++ { + y := uint64(0) for { - xok := false + yok := false for _, ifd := range ifds { - if x < ifd.maxx { - xok = true + if y < ifd.maxy { + yok = true break } } - if !xok { + if !yok { break } - for _, ifd := range ifds { - if x >= ifd.minx && x < ifd.maxx && y >= ifd.miny && y < ifd.maxy { - ch <- tile{ - ifd: ifd, - x: x - ifd.minx, - y: y - ifd.miny, + x := uint64(0) + for { + xok := false + for _, ifd := range ifds { + if x < ifd.maxx { + xok = true + break } } + if !xok { + break + } + for _, ifd := range ifds { + if x >= ifd.minx && x < ifd.maxx && y >= ifd.miny && y < ifd.maxy { + ch <- tile{ + ifd: ifd, + plane: plane, + x: x - ifd.minx, + y: y - ifd.miny, + } + } + } + x++ } - x++ + y++ } - y++ } - for i := 1; i < len(d); i++ { + for i := 0; i < len(d)-1; i++ { ifds = d[i] for _, ifd := range ifds { for y := uint64(0); y < ifd.ntilesy; y++ { for x := uint64(0); x < ifd.ntilesx; x++ { - ch <- tile{ifd: ifd, x: x, y: y} + for p := uint64(0); p < ifd.nplanes; p++ { + ch <- tile{ifd: ifd, x: x, y: y, plane: p} + } } } } diff --git a/testdata/main.go b/testdata/main.go new file mode 100644 index 0000000..af9cb11 --- /dev/null +++ b/testdata/main.go @@ -0,0 +1,73 @@ +package main + +import ( + "fmt" + + "github.com/airbusgeo/godal" +) + +func main() { + godal.RegisterInternalDrivers() + generate1("s1.tif", [4]byte{1, 2, 3, 4}) + generate1("s2.tif", [4]byte{5, 6, 7, 8}) + generate2("p1.tif", "PIXEL", [8]byte{1, 2, 3, 4, 5, 6, 7, 8}) + generate2("p2.tif", "PIXEL", [8]byte{9, 10, 11, 12, 13, 14, 15, 16}) + generate2("b1.tif", "BAND", [8]byte{1, 2, 3, 4, 5, 6, 7, 8}) + generate2("b2.tif", "BAND", [8]byte{9, 10, 11, 12, 13, 14, 15, 16}) +} + +func fillbuf(buf []byte, val byte) { + for i := range buf { + buf[i] = val + } +} +func generate1(fname string, vals [4]byte) { + ds, _ := godal.Create(godal.GTiff, fname, 1, godal.Byte, 128, 128, godal.CreationOption( + "TILED=YES", "BLOCKXSIZE=64", "BLOCKYSIZE=64", + )) + ds.SetGeoTransform([6]float64{0, 0.001, 0, 0, 0, -0.001}) + sr, _ := godal.NewSpatialRefFromEPSG(4326) + ds.SetSpatialRef(sr) + sr.Close() + buf := make([]byte, 64*64) + fillbuf(buf, vals[0]) + ds.Write(0, 0, buf, 64, 64) + fillbuf(buf, vals[1]) + ds.Write(64, 0, buf, 64, 64) + fillbuf(buf, vals[2]) + ds.Write(0, 64, buf, 64, 64) + fillbuf(buf, vals[3]) + ds.Write(64, 64, buf, 64, 64) + ds.BuildOverviews(godal.Levels(2)) + ds.Close() +} +func generate2(fname string, interleave string, vals [8]byte) { + ds, _ := godal.Create(godal.GTiff, fname, 2, godal.Byte, 128, 128, godal.CreationOption( + "TILED=YES", "BLOCKXSIZE=64", "BLOCKYSIZE=64", fmt.Sprintf("INTERLEAVE=%s", interleave), + )) + ds.SetGeoTransform([6]float64{0, 0.001, 0, 0, 0, -0.001}) + sr, _ := godal.NewSpatialRefFromEPSG(4326) + ds.SetSpatialRef(sr) + sr.Close() + buf := make([]byte, 64*64) + bnds := ds.Bands() + fillbuf(buf, vals[0]) + bnds[0].Write(0, 0, buf, 64, 64) + fillbuf(buf, vals[1]) + bnds[0].Write(64, 0, buf, 64, 64) + fillbuf(buf, vals[2]) + bnds[0].Write(0, 64, buf, 64, 64) + fillbuf(buf, vals[3]) + bnds[0].Write(64, 64, buf, 64, 64) + + fillbuf(buf, vals[4]) + bnds[1].Write(0, 0, buf, 64, 64) + fillbuf(buf, vals[5]) + bnds[1].Write(64, 0, buf, 64, 64) + fillbuf(buf, vals[6]) + bnds[1].Write(0, 64, buf, 64, 64) + fillbuf(buf, vals[7]) + bnds[1].Write(64, 64, buf, 64, 64) + ds.BuildOverviews(godal.Levels(2)) + ds.Close() +}