Uneven raster cell size outputs

[pierotofy]

Currently outputs from ODM have slightly uneven cell raster sizes:

>>> import rasterio; rasterio.open("/datasets/brighton2/odm_orthophoto/odm_orthophoto.tif").transform
Affine(0.04999088883357982, 0.0, 576647.5603842067,
       0.0, -0.0499903086638842, 5188220.140716478)
>>> 

This could be improved.

[originlake]

I think the odm_orthophoto renders the orthophoto in the exact raster size provided. It calculates the orthophoto height and width by https://github.com/OpenDroneMap/odm_orthophoto/blob/6a7e81b7dda75e2b71d762f4b6ac89ed97e118c7/src/OdmOrthoPhoto.cpp#L383-L389 Noticed that the ceiling operation there will change the bounds of the orthophoto. But when output corner file, it still uses the original bounds. And later, gdal_translate will inject the old bounds into the orthophoto, which makes the cell size uneven. Therefore, after the ceiling, the bounds should be changed to

bounds.xMax = bounds.xMin + width/resolution_
bounds.yMin = bounds.yMax - height/resolution_

But this still won't fix the problem, because gdal_translate has to calculate the scale from bounds again, and floating-point error is inevitable in such a way. I corrected one orthophoto with the above bounds, it generates the transform like below

Data axis to CRS axis mapping: 1,2
Origin = (728469.337334620999172,5118402.260040946304798)
Pixel Size = (0.050000000000008,-0.050000000000053)

If you can confirm that odm_orthophoto renders the exact raster size as provided, we can directly write the scale value using rasterio (I don't find options to do it with gdal_translate).

[pierotofy]

Thanks for looking into this. Like you say, the core issue is that we manually set the bounds, and precision issues are inevitable.

gdalwarp is able to adjust/warp the cells and create evenly sized cells, but it's slower than gdal_translate.

A really good improvement would be to skip the entire odm_orthophoto --> gdal_translate chain and generate a georeferenced raster directly in odm_orthophoto via calls to GDAL's API.

[originlake]

My point is that odm_orthophoto has already created evenly sized cells, based on how it calculates the orthophoto width and height. So we don't need to adjust/warp the cells, the only problem is then how we properly write the transform data into the orthophoto.

The current approach have two issues

1. The bounds provided to gdal_translate are wrong, it doesn't use bounds after the ceiling operation
2. When using gdal_translate, the scale is calculated with (lrx - ulx)/width, this won't recover the exact resolution provided. The result will be very closed though, 0.050000000000008 vs 0.05 for example.

Issue 1 can be solved by the formula in my last comment Issue 2 can be solved by writing scale and origin directly into the orthophoto, using GDAL API I suppose.

I can prepare a pr for fixing issue 1 firstly if my understanding of odm_orthophto is correct, need your confirmation.

[pierotofy]

I can prepare a pr for fixing issue 1 firstly if my understanding of odm_orthophto is correct, need your confirmation.

Yes I think this is the correct approach (fix 1, then 2).

Careful where you adjust the bounds, I would try to calculate them (with the proposed correction) in computeBoundsForModel or somewhere right after it, in multispectral datasets the model bounds must match (https://github.com/OpenDroneMap/odm_orthophoto/blob/6a7e81b7dda75e2b71d762f4b6ac89ed97e118c7/src/OdmOrthoPhoto.cpp#L365-L380) and you'll get errors if they don't. (Here's an example dataset you can test against to make sure it works): https://github.com/pierotofy/drone_dataset_micasense

[pierotofy]

Forgot to also mention, this issue also affects DSM/DTM outputs (odm_dem/dsm.tif and odm_dem/dtm.tif), I haven't thought about how to fix those yet.. :bulb:

[originlake]

I made some changes to odm_orthophoto, it will calculate the correct bounds, and it can also use GDAL API to write GeoTransfrom and CRS to the geotiff. But then I noticed that gdal_translate will also be used for cutline and feather. Which makes the cell sizes not precise again.

The same thing applied for DSM, the PDAL pipeline generates all the tiff with precise cell sizes, but later GDAL related operations will cause slight precision loss. (It might have something to do with the ceiling and splitting approach, I need to check)

But the precision level is pretty much the float type's level, so I think it should be fine.

[originlake]

https://github.com/OpenDroneMap/odm_orthophoto/pull/1 Updated odm_orthophoto, it will calculate bounds as described above, and work with multispectral data as well. Also added 2 optional arguments to allow directly write geotransform and crs in rendered orthophoto, which will provide better precision, but need to update ODM code. I don't know the plan for staging the odm_orthophoto so I don't modify SuperBuild though.

Before

>>> rio.open('odm_orthophoto.original.tif').transform
Affine(0.009999747968734193, 0.0, 499007.48860168457,
       0.0, -0.009999927561680009, 4857499.880813599)
>>> rio.open('odm_orthophoto.tif').transform
Affine(0.009999747968733269, 0.0, 499009.29855606693,
       0.0, -0.009999927561657001, 4857498.260825333)

Bounds fix

>>> rio.open('odm_orthophoto.original.tif').transform
Affine(0.010000000339581568, 0.0, 499007.48860168457,
       0.0, -0.009999999583488338, 4857499.880813599)
>>> rio.open('odm_orthophoto.tif').transform
Affine(0.010000000339583, 0.0, 499009.29860174604,
       0.0, -0.009999999583505833, 4857498.2608136665)

Directly set GeoTransform

>>> rio.open('odm_orthophoto.original.tif').transform
Affine(0.01, 0.0, 499007.48860168457,
       0.0, -0.01, 4857499.880813599)
>>> rio.open('odm_orthophoto.tif').transform
Affine(0.009999999999999279, 0.0, 499009.29860168457,
       0.0, -0.010000000000008606, 4857498.2608135985)

[pierotofy]

Noticed a (possible?) issue with #1499: even though the cell size is now more uniform, the orthophotos seem to have shifted:

I overlayed them to the DSM, and I think the shift worsens the georeferencing accuracy (the rooftop moves further away from the DSM edge):

I would have expected the overall position of the orthophoto not to shift?

[originlake]

Is it due to scale or shift (is the top left corner of the orthophoto shifted as well)? Seems around 1/4 pixel shift, multiplied by 0.05m is 0.0125m. Assume before was 0.04999088883357982 and later is 0.05, this error could happen when larger than 1371 pixels. If so, I need to have a deeper look at the odm_orthophoto rendering approach.

[originlake]

Yes, I can confirm it's caused by the raster size change. Shifting will increase when the pixel gets far away from the top left corner, no more than 1 pixel. But I'm pretty confident that odm_orthophoto renders texture exactly as given cell size. The simplified conversion from x/y to its corresponding row/column can be represented below (correct me if I'm wrong) and there is no warping around the bounding box.

# given x, y from mesh model
res = 0.05 # cell size in meter
res_ = 1/res # pixels/meter
col = floor((yMax - y) * res_ + 0.5)
row = floor((x - xMin) * res_ + 0.5)

When investigating this, I learned that in GeoTIFF, [0,0] pixel's coordinate has a half pixel shift of the geotransform offsets. When georeferencing the GeoTIFF, we seem not to consider this, I'm not sure if this is relevant.

>>> import rasterio as rio
>>> dst = rio.open('odm_orthophoto.tif')
>>> dst.transform
Affine(0.5, 0.0, 576647.6989822388,
       0.0, -0.5, 5188227.042106628)
>>> dst.xy(0,0)
(576647.9489822388, 5188226.792106628)
>>> dst.xy(-0.5, -0.5)
(576647.6989822388, 5188227.042106628)

[pierotofy]

I think rasterio lets you specify how you want to interpret the xy coordinates (either corner or center, by default it's center), so I don't think it's GeoTIFF-specific thing:

https://rasterio.readthedocs.io/en/latest/api/rasterio.transform.html#rasterio.transform.TransformMethodsMixin.xy

offset (str, optional) – Determines if the returned coordinates are for the center of the pixel or for a corner. Default: center

[originlake]

The problem is then how PDAL defines the origin when it's creating a raster if we are comparing it with the DSM data. And depends on the different implementations, there could be an alignment issue between dem and orthophoto. This might be not related to the raster cell size problem, but if the comparison is between dem and orthophoto, should make it clear. I can tell several differences between pdal and odm_orthophoto.

1. odm_orthophoto aligns the origin at the top left corner while pdal aligns at the bottom left corner
2. odm_orthophoto selects grid center at i,j while PDAL selects grid center at i+0.5, j+0.5source (this seems to cause the difference of the xy interpretation)
3. because of 1, the offsets are calculated differently

These could potentially cause problems and make the above comparison less valid? I'll try to replicate the way PDAL renders raster and see how the results look like

[originlake]

Some updates. I haven't started working on replicating the pdal rendering approach in odm_orthophoto, but managed to render an "RGB" version of DSM. Basically what I did is updating the dsm render pipeline to use RGB value instead of Z value. https://github.com/OpenDroneMap/ODM/blob/5dd0859f4760f3c5a463e5b066dadfde8a53d96a/opendm/dem/pdal.py#L56-L63 Changed them to

d = {
        'type': 'writers.gdal',
        'resolution': resolution,
        'radius': radius,
        'filename': filename,
        'output_type': output_type,
        'data_type': 'uint8',
        'dimension': 'BLUE'
    }

dimensions will be RED, GREEN, and BLUE. Generate "dsm" 3 times with each color. Then use gdal_merge.py to merge bands. And this is what I got, it looks kinda cute : ) By using this to compare with the orthophoto generated from master and 288, I think 288 is closer.

288 vs rgb_dsm

master vs rgb_dsm

I put the generated rgb_dsm in the below link. There are two versions, one is rgb_render_by_odm_dem.tif which is generated by the above approach, and the other version is rgb_render_by_script.tif generated by PDAL directly (there is also a bash script in the folder, that contains the commands used). Both of them are closer to orthophoto generated in branch 288.

https://drive.google.com/drive/folders/1aVE5gxOvSktPxOeOizZcsDCJj52bmM34?usp=sharing

[originlake]

2. odm_orthophoto selects grid center at i,j while PDAL selects grid center at i+0.5, j+0.5source (this seems to cause the difference of the xy interpretation)

I was wrong about this, odm_orthophoto uses i,j to obtain triangle boundary but, when actually render pixel color, it uses i+0.5, j+0.5

getBarycentricCoordinates(v1, v2, v3, static_cast<float>(cq)+0.5f, static_cast<float>(rq)+0.5f, l1, l2, l3);

Therefore, odm_orthophoto also uses pixel center as xy. Same as rasterio defaults or pdal. There is no difference for the corner or center assumptions.

I believe the use PDAL directly rendered RGB "DSM" to compare is more straightforward than DSM. The reason why DSM seems to have offsets might be related to the data_type=max and radiuses. I could imagine some edge dilation happens for objects that have higher elevations with this setting. I can create similar DSM "offsets" by reducing the radius in odm_dem by half.

[pierotofy]

Thanks for looking further into this!

Yes PDAL will always need a radius parameter to render results using writers.gdal, so whether it's RGB or elevation, tweaking the radius parameter will alter results.

I suppose a better source for ground truth should be gdalwarp. For example, I noticed I can generate uniform raster cell sizes with a simple call to gdalwarp:

gdalwarp odm_orthophoto.tif warped.tif

>>> import rasterio
>>> rasterio.open('odm_orthophoto.tif').transform
Affine(0.04999806992055534, 0.0, 576649.6427728771,
       0.0, -0.049992216122986016, 5188210.877415851)
>>> rasterio.open('warped.tif').transform
Affine(0.04999575313770262, 0.0, 576649.6427728771,
       0.0, -0.04999575313770262, 5188210.877415851)

When opened in QGIS, the two files are slightly different, but the pixel shift is much less apparent.

[originlake]

I still believe the uneven problem is caused by the problem I mentioned in the first comment. Given an example pointcloud corners

xMin = -6.78768844604492188e+01 
yMin = -6.51020584106445312e+01 
xMax = 7.20856628417968750e+01
yMax = 7.20272445678710938e+01
resolution = 0.05

Odm_orthophoto firstly calculates the image size by

height = ceil((yMax-yMin)/resolution) # 2743
width = ceil((xMax-xMin)/resolution) # 2800

Then it renders pixel size exactly at 0.05m/pixel. Can be proved by that it transforms the points by the below affine matrix so that the points are directly located inside the image grid.

1/0.05 , 0 , -xMin 0 , -1/0.05 , yMax 0 , 0 , 1

The ceil operation will extend the bounds, but if we still use the old bounds to wrap the orthophoto, it will cause errors. And because x and y errors are different, which makes the cell size uneven. The errors are not big, they are limited by the resolution, but half of the resolution error means half a pixel offset at the bottom right corner of the image.

diff_y = height * resolution - (yMax-yMin) # 0.020697021484380684
diff_x = width * resolution - (xMax-xMin) # 0.03745269775390625
# calculate cell size by gdal_translate, the results match the cell size in odm_orthophoto.original.tif
res_y = (yMin - yMax)/height # -0.049992454603906535
res_x = (xMax - xMin)/width # 0.049986624036516464

If the image is rendered with the extended bounds(larger) but wrapped with the original bounds(smaller). This will shrink the cell size a little bit. The difference for each cell is tiny, but it will accumulate from tl corner to br corner, this is the reason for the offsets. And I do think the approach I use PDAL to generate orthophoto provides valid proof that my fix is correct.

Anyway, this is all I got from the recent investigation, I hope this can persuade you.

[pierotofy]

I started looking at this problem a bit from the DEM generation.

The cells become uneven during the crop stage (here: https://github.com/OpenDroneMap/ODM/blob/master/opendm/cropper.py#L53), I suspect because the crop polygon is not aligned to properly round boundaries. Interestingly one can pass -tr <xres> <yres> and get even raster cells, but a 1/2 pixel shift happens.

Maybe the way to do this would be to simply add a new --even-rastersize flag explaining that the pixel shift will happen as a tradeoff for having exact raster cell sizes.

[originlake]

Sounds good!