diff options
| -rw-r--r-- | data/Makefile.am | 7 | ||||
| -rw-r--r-- | data/tests/tinshift_crs_implicit.json | 46 | ||||
| -rw-r--r-- | data/tests/tinshift_simplified_kkj_etrs.json | 50 | ||||
| -rw-r--r-- | data/tests/tinshift_simplified_n60_n2000.json | 50 | ||||
| -rw-r--r-- | data/triangulation.schema.json | 169 | ||||
| -rw-r--r-- | docs/Makefile | 1 | ||||
| -rw-r--r-- | docs/source/operations/transformations/index.rst | 1 | ||||
| -rw-r--r-- | docs/source/operations/transformations/tinshift.rst | 201 | ||||
| -rwxr-xr-x | scripts/reformat_cpp.sh | 5 | ||||
| -rw-r--r-- | src/Makefile.am | 5 | ||||
| -rw-r--r-- | src/lib_proj.cmake | 1 | ||||
| -rw-r--r-- | src/pj_list.h | 1 | ||||
| -rw-r--r-- | src/quadtree.hpp | 257 | ||||
| -rw-r--r-- | src/transformations/tinshift.cpp | 133 | ||||
| -rw-r--r-- | src/transformations/tinshift.hpp | 257 | ||||
| -rw-r--r-- | src/transformations/tinshift_exceptions.hpp | 52 | ||||
| -rw-r--r-- | src/transformations/tinshift_impl.hpp | 555 | ||||
| -rw-r--r-- | test/CMakeLists.txt | 1 | ||||
| -rw-r--r-- | test/gie/Makefile.am | 8 | ||||
| -rw-r--r-- | test/gie/tinshift.gie | 50 | ||||
| -rw-r--r-- | test/unit/CMakeLists.txt | 11 | ||||
| -rw-r--r-- | test/unit/Makefile.am | 9 | ||||
| -rw-r--r-- | test/unit/test_tinshift.cpp | 211 |
23 files changed, 2077 insertions, 4 deletions
diff --git a/data/Makefile.am b/data/Makefile.am index b48f53fa..bdfabe0f 100644 --- a/data/Makefile.am +++ b/data/Makefile.am @@ -4,7 +4,8 @@ pkgdata_DATA = proj.ini GL27 nad.lst nad27 nad83 world other.extra \ CH \ ITRF2000 ITRF2008 ITRF2014 proj.db \ projjson.schema.json \ - deformation_model.schema.json + deformation_model.schema.json \ + triangulation.schema.json SQL_ORDERED_LIST = sql/begin.sql \ sql/proj_db_table_defs.sql \ @@ -45,6 +46,7 @@ EXTRA_DIST = proj.ini GL27 nad.lst nad27 nad83 \ ITRF2000 ITRF2008 ITRF2014 \ projjson.schema.json \ deformation_model.schema.json \ + triangulation.schema.json \ CMakeLists.txt \ tests/test_nodata.gtx \ tests/test_vgrid_bigendian_bigtiff.tif \ @@ -116,6 +118,9 @@ EXTRA_DIST = proj.ini GL27 nad.lst nad27 nad83 \ tests/simple_model_wrap_west.json \ tests/simple_model_wrap_west.tif \ tests/simple_model_projected.json \ + tests/tinshift_crs_implicit.json \ + tests/tinshift_simplified_kkj_etrs.json \ + tests/tinshift_simplified_n60_n2000.json \ generate_all_sql_in.cmake sql_filelist.cmake \ $(SQL_ORDERED_LIST) diff --git a/data/tests/tinshift_crs_implicit.json b/data/tests/tinshift_crs_implicit.json new file mode 100644 index 00000000..e1939fd0 --- /dev/null +++ b/data/tests/tinshift_crs_implicit.json @@ -0,0 +1,46 @@ +{ + "file_type": "triangulation_file", + "format_version": "1.0", + "name": "Name", + "version": "Version", + "publication_date": "2018-07-01T00:00:00Z", + "license": "Creative Commons Attribution 4.0 International", + "description": "Test triangulation", + "authority": { + "name": "Authority name", + "url": "http://example.com", + "address": "Adress", + "email": "test@example.com" + }, + "links": [ + { + "href": "https://example.com/about.html", + "rel": "about", + "type": "text/html", + "title": "About" + }, + { + "href": "https://example.com/download", + "rel": "source", + "type": "application/zip", + "title": "Authoritative source" + }, + { + "href": "https://creativecommons.org/licenses/by/4.0/", + "rel": "license", + "type": "text/html", + "title": "Creative Commons Attribution 4.0 International license" + }, + { + "href": "https://example.com/metadata.xml", + "rel": "metadata", + "type": "application/xml", + "title": " ISO 19115 XML encoded metadata regarding the deformation model" + } + ], + "transformed_components": [ "horizontal" ], + "vertices_columns": [ "source_x", "source_y", "target_x", "target_y" ], + "triangles_columns": [ "idx_vertex1", "idx_vertex2", "idx_vertex3" ], + "vertices": [ [2,49,2.1,49.1], [3,50,3.1,50.1], [2, 50, 2.1,50.1] ], + "triangles": [ [0, 1, 2] ] +} diff --git a/data/tests/tinshift_simplified_kkj_etrs.json b/data/tests/tinshift_simplified_kkj_etrs.json new file mode 100644 index 00000000..7d46107f --- /dev/null +++ b/data/tests/tinshift_simplified_kkj_etrs.json @@ -0,0 +1,50 @@ +{ + "file_type": "triangulation_file", + "format_version": "1.0", + "name": "Name", + "version": "Version", + "publication_date": "2018-07-01T00:00:00Z", + "license": "Creative Commons Attribution 4.0 International", + "description": "Test triangulation", + "authority": { + "name": "Authority name", + "url": "http://example.com", + "address": "Adress", + "email": "test@example.com" + }, + "links": [ + { + "href": "https://example.com/about.html", + "rel": "about", + "type": "text/html", + "title": "About" + }, + { + "href": "https://example.com/download", + "rel": "source", + "type": "application/zip", + "title": "Authoritative source" + }, + { + "href": "https://creativecommons.org/licenses/by/4.0/", + "rel": "license", + "type": "text/html", + "title": "Creative Commons Attribution 4.0 International license" + }, + { + "href": "https://example.com/metadata.xml", + "rel": "metadata", + "type": "application/xml", + "title": " ISO 19115 XML encoded metadata regarding the triangulation" + } + ], + "input_crs": "EPSG:2393", + "output_crs": "EPSG:3067", + "transformed_components": [ "horizontal" ], + "vertices_columns": [ "source_x", "source_y", "target_x", "target_y" ], + "triangles_columns": [ "idx_vertex1", "idx_vertex2", "idx_vertex3" ], + "vertices": [ [3244102.707, 6693710.937, 244037.137, 6690900.686], + [3205290.722, 6715311.822, 205240.895, 6712492.577], + [3218328.492, 6649538.429, 218273.648, 6646745.973] ], + "triangles": [ [0, 1, 2] ] +} diff --git a/data/tests/tinshift_simplified_n60_n2000.json b/data/tests/tinshift_simplified_n60_n2000.json new file mode 100644 index 00000000..f6d3d4aa --- /dev/null +++ b/data/tests/tinshift_simplified_n60_n2000.json @@ -0,0 +1,50 @@ +{ + "file_type": "triangulation_file", + "format_version": "1.0", + "name": "Name", + "version": "Version", + "publication_date": "2018-07-01T00:00:00Z", + "license": "Creative Commons Attribution 4.0 International", + "description": "Test triangulation", + "authority": { + "name": "Authority name", + "url": "http://example.com", + "address": "Adress", + "email": "test@example.com" + }, + "links": [ + { + "href": "https://example.com/about.html", + "rel": "about", + "type": "text/html", + "title": "About" + }, + { + "href": "https://example.com/download", + "rel": "source", + "type": "application/zip", + "title": "Authoritative source" + }, + { + "href": "https://creativecommons.org/licenses/by/4.0/", + "rel": "license", + "type": "text/html", + "title": "Creative Commons Attribution 4.0 International license" + }, + { + "href": "https://example.com/metadata.xml", + "rel": "metadata", + "type": "application/xml", + "title": " ISO 19115 XML encoded metadata regarding the tirangulation" + } + ], + "input_crs": "EPSG:2393+5717", + "output_crs": "EPSG:2393+5941", + "transformed_components": [ "vertical" ], + "vertices_columns": [ "source_x", "source_y", "source_z", "target_z" ], + "triangles_columns": [ "idx_vertex1", "idx_vertex2", "idx_vertex3" ], + "vertices": [ [3188607.0, 6688748.0, 23.123, 23.4133], + [3184981.0, 6725255.0, 8.044, 8.34499], + [3220912.0, 6699508.0, 1.724, 2.0101] ], + "triangles": [ [0, 1, 2] ] +} diff --git a/data/triangulation.schema.json b/data/triangulation.schema.json new file mode 100644 index 00000000..ee5cae0c --- /dev/null +++ b/data/triangulation.schema.json @@ -0,0 +1,169 @@ +{ + "$schema": "http://json-schema.org/draft-07/schema#", + "description": "Schema for triangulation based transformation", + "type": "object", + "properties": { + "file_type": { + "type": "string", + "enum": [ + "triangulation_file" + ], + "description": "File type. Always \"triangulation_file\"" + }, + "format_version": { + "type": "string", + "enum": [ + "1.0" + ] + }, + "name": { + "type": "string", + "description": "A brief descriptive name of the triangulation" + }, + "version": { + "type": "string", + "description": "A string identifying the version of the triangulation. The format for specifying version will be defined by the agency responsible for the triangulation" + }, + "publication_date": { + "$ref": "#/definitions/datetime", + "description": "The date on which this version of the triangulation was published (or possibly the date on which it takes effect?)" + }, + "license": { + "type": "string", + "description": "License under which the file is published" + }, + "description": { + "type": "string", + "description": "A text description of the file" + }, + "authority": { + "type": "object", + "description": "Basic information about the agency responsible for the data set", + "properties": { + "name": { + "type": "string", + "description": "The name of the agency" + }, + "url": { + "type": "string", + "description": "The url of the agency website", + "format": "uri" + }, + "address": { + "type": "string", + "description": "The postal address of the agency" + }, + "email": { + "type": "string", + "description": "An email contact address for the agency", + "format": "email" + } + }, + "required": [ + "name" + ], + "additionalProperties": false + }, + "links": { + "type": "array", + "description": "Links to related information", + "items": { + "type": "object", + "properties": { + "href": { + "type": "string", + "description": "The URL holding the information", + "format": "uri" + }, + "rel": { + "type": "string", + "description": "The relationship to the dataset. Proposed relationships are:\n- \"about\": a web page for human consumption describing the model\n- \"source\": the authoritative source data from which the triangulation is built.\n- \"metadata\": ISO 19115 XML metadata regarding the triangulation." + }, + "type": { + "type": "string", + "description": "MIME type" + }, + "title": { + "type": "string", + "description": "Description of the link" + } + }, + "required": [ + "href" + ], + "additionalProperties": false + } + }, + "input_crs": { + "$ref": "#/definitions/crs", + "description": "String identifying the CRS of source coordinates in the vertices. Typically \"EPSG:XXXX\". If the transformation is for vertical component, this should be the code for a compound CRS (can be EPSG:XXXX+YYYY where XXXX is the code of the horizontal CRS and YYYY the code of the vertical CRS). For example, for the KKJ->ETRS89 transformation, this is EPSG:2393 (\"KKJ / Finland Uniform Coordinate System\"). The input coordinates are assumed to be passed in the \"normalized for visualisation\" / \"GIS friendly\" order, that is longitude, latitude for geographic coordinates and easting, northing for projected coordinates." + }, + "output_crs": { + "$ref": "#/definitions/crs", + "description": "String identifying the CRS of target coordinates in the vertices. Typically \"EPSG:XXXX\". If the transformation is for vertical component, this should be the code for a compound CRS (can be EPSG:XXXX+YYYY where XXXX is the code of the horizontal CRS and YYYY the code of the vertical CRS). For example, for the KKJ->ETRS89 transformation, this is EPSG:3067 (\"ETRS89 / TM35FIN(E,N)\"). The output coordinates will be returned in the \"normalized for visualisation\" / \"GIS friendly\" order, that is easting, that is longitude, latitude for geographic coordinates and easting, northing for projected coordinates." + }, + "transformed_components": { + "type": "array", + "description": "Specify which component of the coordinates are transformed. Either \"horizontal\", \"vertical\" or both", + "minItems": 1, + "maxItems": 2, + "items": { + "type": "string", + "enum": [ + "horizontal", + "vertical" + ] + } + }, + "vertices_columns": { + "type": "array", + "description": "Specify the name of the columns of the rows in the \"vertices\" array. There must be exactly as many elements in \"vertices_columns\" as in a row of \"vertices\". The following names have a special meaning: \"source_x\", \"source_y\", \"target_x\", \"target_y\", \"source_z\", \"target_z\" and \"offset_z\". \"source_x\" and \"source_y\" are compulsory. \"source_x\" is for the source longitude (in degree) or easting. \"source_y\" is for the source latitude (in degree) or northing. \"target_x\" and \"target_y\" are compulsory when \"horizontal\" is specified in \"transformed_components\". (\"source_z\" and \"target_z\") or \"offset_z\" are compulsory when \"vertical\" is specified in \"transformed_components\".", + "minItems": 3, + "items": { + "type": "string" + } + }, + "triangles_columns": { + "type": "array", + "description": "Specify the name of the columns of the rows in the \"triangles\" array. There must be exactly as many elements in \"triangles_columns\" as in a row of \"triangles\". The following names have a special meaning: \"idx_vertex1\", \"idx_vertex2\", \"idx_vertex3\". They are compulsory.", + "minItems": 3, + "items": { + "type": "string" + } + }, + "vertices": { + "type": "array", + "description": "an array whose items are themselves arrays with as many columns as described in \"vertices_columns\"", + "items": { + "type": "array" + } + }, + "triangles": { + "type": "array", + "description": "an array whose items are themselves arrays with as many columns as described in \"triangles_columns\". The value of the \"idx_vertexN\" columns must be indices (between 0 and len(\"vertices\"-1) of items of the \"vertices\" array", + "items": { + "type": "array" + } + } + }, + "required": [ + "file_type", + "format_version", + "transformed_components", + "vertices_columns", + "triangles_columns", + "vertices", + "triangles" + ], + "additionalProperties": false, + "definitions": { + "crs": { + "type": "string" + }, + "datetime": { + "type": "string", + "format": "date-time", + "pattern": "^[0-9]{4}-[0-9]{2}-[0-9]{2}T[0-9]{2}:[0-9]{2}:[0-9]{2}Z$" + } + } +} diff --git a/docs/Makefile b/docs/Makefile index c0801a9b..a63498ab 100644 --- a/docs/Makefile +++ b/docs/Makefile @@ -66,6 +66,7 @@ html: .doxygen_up_to_date @sed "s/<em class=\"property\">namespace <\/em>//g" < $(BUILDDIR)/html/development/reference/cpp/cpp_general.html > $(BUILDDIR)/html/development/reference/cpp/cpp_general.html.tmp @mv $(BUILDDIR)/html/development/reference/cpp/cpp_general.html.tmp $(BUILDDIR)/html/development/reference/cpp/cpp_general.html @cp -r ../schemas $(BUILDDIR)/html/schemas + @cp ../data/triangulation.schema.json $(BUILDDIR)/html/schemas @echo @echo "Build finished. The HTML pages are in $(BUILDDIR)/html." diff --git a/docs/source/operations/transformations/index.rst b/docs/source/operations/transformations/index.rst index 9f7900f1..a413b8f6 100644 --- a/docs/source/operations/transformations/index.rst +++ b/docs/source/operations/transformations/index.rst @@ -19,5 +19,6 @@ systems are based on different datums. molodensky molobadekas hgridshift + tinshift vgridshift xyzgridshift diff --git a/docs/source/operations/transformations/tinshift.rst b/docs/source/operations/transformations/tinshift.rst new file mode 100644 index 00000000..c29d233a --- /dev/null +++ b/docs/source/operations/transformations/tinshift.rst @@ -0,0 +1,201 @@ +.. _tinshift: + +================================================================================ +Triangulated Irregular Network based transformation +================================================================================ + +.. versionadded:: 7.2.0 + ++---------------------+--------------------------------------------------------------------+ +| **Alias** | tinshift | ++---------------------+--------------------------------------------------------------------+ +| **Input type** | Projected or geographic coordinates (horizontal), meters (vertical)| ++---------------------+--------------------------------------------------------------------+ +| **Output type** | Projected or geographic coordinates (horizontal), meters (vertical)| ++---------------------+--------------------------------------------------------------------+ +| **Domain** | 2D or 3D | ++---------------------+--------------------------------------------------------------------+ +| **Available forms** | Forward and inverse | ++---------------------+--------------------------------------------------------------------+ + + +The ``tinshift`` transformation takes one mandatory +argument, ``file``, that points to a JSON file, which contains the +triangulation and associated metadata. Input and output coordinates must be +geographic or projected coordinates. +Depending on the content of the JSON file, horizontal, vertical or both +components of the coordinates may be transformed. + +The transformation is invertible, with the same computational complexity than +the forward transformation. + +Parameters +------------------------------------------------------------------------------- + +Required ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ + +.. option:: +file=<filename> + + Filename to the JSON file for the TIN. + + +Example +------------------------------------------------------------------------------- + +Transformating a point with the Finland EPSG:2393 ("KKJ / Finland Uniform +Coordinate System") projected CRS to EPSG:3067 ("ETRS89 / TM35FIN(E,N)") + +:: + + $ echo 3210000.0000 6700000.0000 0 2020 | cct +proj=tinshift +file=./triangulation_kkj.json + + 209948.3217 6697187.0009 0.0000 2020 + + +Algorithm ++++++++++ + +Internally, ``tinshift`` ingest the whole file into memory. It is considered that +triangulation should be small enough for that. + +When a point is transformed, one must find the triangle into which it falls into. +Instead of iterating over all triangles, we build a in-memory quadtree to speed-up +the identification of candidates triangles. + +To determine if a point falls into a triangle, one computes its 3 +`barycentric coordinates <https://en.wikipedia.org/wiki/Barycentric_coordinate_system#Conversion_between_barycentric_and_Cartesian_coordinates>`_ +from its projected coordinates, :math:`\lambda_i` for :math:`i=1,2,3`. +They are real values (in the [0,1] range for a point inside the triangle), +giving the weight of each of the 3 vertices of the triangles. + +Once those weights are known, interpolating the target horizontal +coordinate is a matter of doing the linear combination of those weights with +the target horizontal coordinates at the 3 vertices of the triangle (:math:`Xt_i` and :math:`Yt_i`): + +.. math:: + + X_{target} = Xt_1 * \lambda_1 + Xt_2 * \lambda_2 + Xt_3 * \lambda_3 + + Y_{target} = Yt_1 * \lambda_1 + Yt_2 * \lambda_2 + Yt_3 * \lambda_3 + +This interpolation is exact at the vertices of the triangulation, and has linear properties +inside each triangle. It is completely equivalent to other formulations of +triangular interpolation, such as + +.. math:: + + X_{target} = A + X_{source} * B + Y_{source} * C + + Y_{target} = D + X_{source} * E + Y_{source} * F + +where the A, B, C, D, E, F constants (for a given triangle) are found by solving +the 2 systems of 3 linear equations, constraint by the source and target coordinate pairs +of the 3 vertices of the triangle: + +.. math:: + + Xt_i = A + Xs_i * B + Ys_i * C + + Yt_i = D + Xs_i * E + Ys_i * F + +Similarly for a vertical coordinate transformation, where :math:`Zoff_i` is the vertical +offset at each vertex of the triangle: + +.. math:: + + Z_{target} = Z_{source} + Zoff_1 * \lambda_1 + Zoff_2 * \lambda_2 + Zoff_3 * \lambda_3 + +Constraints on the triangulation +++++++++++++++++++++++++++++++++ + +No check is done on the consistence of the triangulation. It is highly +recommended that triangles do not overlap each other (when considering the +source coordinates or the forward transformation, or the target coordinates for +the inverse transformation), otherwise which triangle will be selected is +unspecified. Besides that, the triangulation does not need to have particular +properties (like being a Delaunay triangulation) + +File format ++++++++++++ + +The triangulation is stored in a text-based format, using JSON as a serialization. + +Below a minimal example, from the KKJ to ETRS89 transformation, with just a +single triangle: + +.. literalinclude:: ../../../../data/tests/tinshift_crs_implicit.json + :language: json + + +So after the generic metadata, we define the input and output CRS (informative +only), and that the transformation affects horizontal components of +coordinates. We name the columns of the ``vertices`` and ``triangles`` arrays. +We defined the source and target coordinates of each vertex, and define a +triangle by referring to the index of its vertices in the ``vertices`` array. + +More formally, the specific items for the triangulation file are: + +input_crs + String identifying the CRS of source coordinates + in the vertices. Typically ``EPSG:XXXX``. If the transformation is for vertical + component, this should be the code for a compound CRS (can be EPSG:XXXX+YYYY + where XXXX is the code of the horizontal CRS and YYYY the code of the vertical CRS). + For example, for the KKJ->ETRS89 transformation, this is EPSG:2393 + (``KKJ / Finland Uniform Coordinate System``). The input coordinates are assumed + to be passed in the "normalized for visualisation" / "GIS friendly" order, + that is longitude, latitude for geographic coordinates and + easting, northing for projected coordinates. + + +output_crs + String identifying the CRS of target coordinates in the vertices. + Typically ``EPSG:XXXX``. If the transformation is for vertical component, + this should be the code for a compound CRS (can be EPSG:XXXX+YYYY where + XXXX is the code of the horizontal CRS and YYYY the code of the vertical CRS). + For example, for the KKJ->ETRS89 transformation, this is EPSG:3067 + (\"ETRS89 / TM35FIN(E,N)\"). The output coordinates will be returned in + the "normalized for visualisation" / "GIS friendly" order, + that is longitude, latitude for geographic coordinates and + easting, northing for projected coordinates. + + +transformed_components + Array which may contain one or two strings: "horizontal" when horizontal + components of the coordinates are transformed and/or "vertical" when the + vertical component is transformed. + + +vertices_columns + Specify the name of the columns of the rows in the ``vertices`` + array. There must be exactly as many elements in ``vertices_columns`` as in a + row of ``vertices``. The following names have a special meaning: ``source_x``, + ``source_y``, ``target_x``, ``target_y``, ``source_z``, ``target_z`` and + ``offset_z``. ``source_x`` and ``source_y`` are compulsory. + ``source_x`` is for the source longitude (in degree) or easting. + ``source_y`` is for the source latitude (in degree) or northing. + ``target_x`` and ``target_y`` are compulsory when ``horizontal`` is specified + in ``transformed_components``. (``source_z`` and ``target_z``) or + ``offset_z`` are compulsory when ``vertical`` is specified in ``transformed_components`` + + +triangles_columns + Specify the name of the columns of the rows in the + ``triangles`` array. There must be exactly as many elements in ``triangles_columns`` + as in a row of ``triangles``. The following names have a special meaning: + ``idx_vertex1``, ``idx_vertex2``, ``idx_vertex3``. They are compulsory. + + +vertices + An array whose items are themselves arrays with as many columns as + described in ``vertices_columns``. + + +triangles + An array whose items are themselves arrays with as many columns as + described in ``triangles_columns``. + The value of the ``idx_vertexN`` columns must be indices + (between 0 and len(``vertices``-1) of items of the ``vertices`` array. + +A `JSON schema <https://proj.org/schemas/triangulation.schema.json>`_ is available +for this file format. diff --git a/scripts/reformat_cpp.sh b/scripts/reformat_cpp.sh index 0c57daba..50a572a1 100755 --- a/scripts/reformat_cpp.sh +++ b/scripts/reformat_cpp.sh @@ -27,6 +27,11 @@ for i in "$TOPDIR"/include/proj/*.hpp "$TOPDIR"/include/proj/internal/*.hpp \ "$TOPDIR"/src/transformations/defmodel_exceptions.hpp \ "$TOPDIR"/src/transformations/defmodel_impl.hpp \ "$TOPDIR"/src/transformations/defmodel.cpp \ + "$TOPDIR"/src/transformations/tinshift.hpp \ + "$TOPDIR"/src/transformations/tinshift_exceptions.hpp \ + "$TOPDIR"/src/transformations/tinshift_impl.hpp \ + "$TOPDIR"/src/transformations/tinshift.cpp \ + "$TOPDIR"/src/quadtree.hpp \ ; do if ! echo "$i" | grep -q "lru_cache.hpp"; then "$SCRIPT_DIR"/reformat.sh "$i"; diff --git a/src/Makefile.am b/src/Makefile.am index 86444df3..68aeadf2 100644 --- a/src/Makefile.am +++ b/src/Makefile.am @@ -201,6 +201,10 @@ libproj_la_SOURCES = \ transformations/defmodel.hpp \ transformations/defmodel_exceptions.hpp \ transformations/defmodel_impl.hpp \ + transformations/tinshift.cpp \ + transformations/tinshift.hpp \ + transformations/tinshift_exceptions.hpp \ + transformations/tinshift_impl.hpp \ \ aasincos.cpp adjlon.cpp \ dmstor.cpp auth.cpp \ @@ -213,6 +217,7 @@ libproj_la_SOURCES = \ release.cpp gauss.cpp \ fileapi.cpp \ generic_inverse.cpp \ + quadtree.hpp \ \ datums.cpp datum_set.cpp transform.cpp \ utils.cpp \ diff --git a/src/lib_proj.cmake b/src/lib_proj.cmake index 3d41fe9a..1cd7f421 100644 --- a/src/lib_proj.cmake +++ b/src/lib_proj.cmake @@ -186,6 +186,7 @@ set(SRC_LIBPROJ_TRANSFORMATIONS transformations/vgridshift.cpp transformations/xyzgridshift.cpp transformations/defmodel.cpp + transformations/tinshift.cpp ) set(SRC_LIBPROJ_ISO19111 diff --git a/src/pj_list.h b/src/pj_list.h index c7e6d209..5b91af9b 100644 --- a/src/pj_list.h +++ b/src/pj_list.h @@ -150,6 +150,7 @@ PROJ_HEAD(gstmerc, "Gauss-Schreiber Transverse Mercator (aka Gauss-Laborde Reuni PROJ_HEAD(tcc, "Transverse Central Cylindrical") PROJ_HEAD(tcea, "Transverse Cylindrical Equal Area") PROJ_HEAD(times, "Times Projection") +PROJ_HEAD(tinshift, "Triangulation based transformation") PROJ_HEAD(tissot, "Tissot Conic") PROJ_HEAD(tmerc, "Transverse Mercator") PROJ_HEAD(tobmerc, "Tobler-Mercator") diff --git a/src/quadtree.hpp b/src/quadtree.hpp new file mode 100644 index 00000000..b744c83d --- /dev/null +++ b/src/quadtree.hpp @@ -0,0 +1,257 @@ +/****************************************************************************** + * Project: PROJ + * Purpose: Implementation of quadtree building and searching functions. + * Derived from shapelib, mapserver and GDAL implementations + * Author: Even Rouault, <even.rouault at spatialys.com> + * + ****************************************************************************** + * Copyright (c) 1999-2008, Frank Warmerdam + * Copyright (c) 2008-2020, Even Rouault <even dot rouault at spatialys.com> + * + * Permission is hereby granted, free of charge, to any person obtaining a + * copy of this software and associated documentation files (the "Software"), + * to deal in the Software without restriction, including without limitation + * the rights to use, copy, modify, merge, publish, distribute, sublicense, + * and/or sell copies of the Software, and to permit persons to whom the + * Software is furnished to do so, subject to the following conditions: + * + * The above copyright notice and this permission notice shall be included + * in all copies or substantial portions of the Software. + * + * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS + * OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, + * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL + * THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER + * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING + * FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER + * DEALINGS IN THE SOFTWARE. + *****************************************************************************/ + +#ifndef QUADTREE_HPP +#define QUADTREE_HPP + +#include "proj/util.hpp" + +#include <functional> +#include <vector> + +//! @cond Doxygen_Suppress + +NS_PROJ_START + +namespace QuadTree { + +/* -------------------------------------------------------------------- */ +/* If the following is 0.5, psNodes will be split in half. If it */ +/* is 0.6 then each apSubNode will contain 60% of the parent */ +/* psNode, with 20% representing overlap. This can be help to */ +/* prevent small objects on a boundary from shifting too high */ +/* up the hQuadTree. */ +/* -------------------------------------------------------------------- */ +constexpr double DEFAULT_SPLIT_RATIO = 0.55; + +/** Describe a rectangle */ +struct RectObj { + double minx = 0; /**< Minimum x */ + double miny = 0; /**< Minimum y */ + double maxx = 0; /**< Maximum x */ + double maxy = 0; /**< Maximum y */ + + /* Returns whether this rectangle is contained by other */ + inline bool isContainedBy(const RectObj &other) const { + return minx >= other.minx && maxx <= other.maxx && miny >= other.miny && + maxy <= other.maxy; + } + + /* Returns whether this rectangles overlaps other */ + inline bool overlaps(const RectObj &other) const { + return minx <= other.maxx && maxx >= other.minx && miny <= other.maxy && + maxy >= other.miny; + } + + /* Returns whether this rectangles contains the specified point */ + inline bool contains(double x, double y) const { + return minx <= x && maxx >= x && miny <= y && maxy >= y; + } + + /* Return whether this rectangles is different from other */ + inline bool operator!=(const RectObj &other) const { + return minx != other.minx || miny != other.miny || maxx != other.maxx || + maxy != other.maxy; + } +}; + +/** Quadtree */ +template <class Feature> class QuadTree { + + struct Node { + RectObj rect{}; /* area covered by this psNode */ + + /* list of shapes stored at this node. */ + std::vector<std::pair<Feature, RectObj>> features{}; + + std::vector<Node> subnodes{}; + + explicit Node(const RectObj &rectIn) : rect(rectIn) {} + }; + + Node root{}; + unsigned nBucketCapacity = 8; + double dfSplitRatio = DEFAULT_SPLIT_RATIO; + + public: + /** Construct a new quadtree with the global bounds of all objects to be + * inserted */ + explicit QuadTree(const RectObj &globalBounds) : root(globalBounds) {} + + /** Add a new feature, with its bounds specified in featureBounds */ + void insert(const Feature &feature, const RectObj &featureBounds) { + insert(root, feature, featureBounds); + } + +#ifdef UNUSED + /** Retrieve all features whose bounds intersects aoiRect */ + void + search(const RectObj &aoiRect, + std::vector<std::reference_wrapper<const Feature>> &features) const { + search(root, aoiRect, features); + } +#endif + + /** Retrieve all features whose bounds contains (x,y) */ + void search(double x, double y, std::vector<Feature> &features) const { + search(root, x, y, features); + } + + private: + void splitBounds(const RectObj &in, RectObj &out1, RectObj &out2) { + // The output bounds will be very similar to the input bounds, + // so just copy over to start. + out1 = in; + out2 = in; + + // Split in X direction. + if ((in.maxx - in.minx) > (in.maxy - in.miny)) { + const double range = in.maxx - in.minx; + + out1.maxx = in.minx + range * dfSplitRatio; + out2.minx = in.maxx - range * dfSplitRatio; + } + + // Otherwise split in Y direction. + else { + const double range = in.maxy - in.miny; + + out1.maxy = in.miny + range * dfSplitRatio; + out2.miny = in.maxy - range * dfSplitRatio; + } + } + + void insert(Node &node, const Feature &feature, + const RectObj &featureBounds) { + if (node.subnodes.empty()) { + // If we have reached the max bucket capacity, try to insert + // in a subnode if possible. + if (node.features.size() >= nBucketCapacity) { + RectObj half1; + RectObj half2; + RectObj quad1; + RectObj quad2; + RectObj quad3; + RectObj quad4; + + splitBounds(node.rect, half1, half2); + splitBounds(half1, quad1, quad2); + splitBounds(half2, quad3, quad4); + + if (node.rect != quad1 && node.rect != quad2 && + node.rect != quad3 && node.rect != quad4 && + (featureBounds.isContainedBy(quad1) || + featureBounds.isContainedBy(quad2) || + featureBounds.isContainedBy(quad3) || + featureBounds.isContainedBy(quad4))) { + node.subnodes.reserve(4); + node.subnodes.emplace_back(Node(quad1)); + node.subnodes.emplace_back(Node(quad2)); + node.subnodes.emplace_back(Node(quad3)); + node.subnodes.emplace_back(Node(quad4)); + + auto features = std::move(node.features); + node.features.clear(); + for (auto &pair : features) { + insert(node, pair.first, pair.second); + } + + /* recurse back on this psNode now that it has apSubNodes */ + insert(node, feature, featureBounds); + return; + } + } + } else { + // If we have sub nodes, then consider whether this object will + // fit in them. + for (auto &subnode : node.subnodes) { + if (featureBounds.isContainedBy(subnode.rect)) { + insert(subnode, feature, featureBounds); + return; + } + } + } + + // If none of that worked, just add it to this nodes list. + node.features.push_back( + std::pair<Feature, RectObj>(feature, featureBounds)); + } + +#ifdef UNUSED + void + search(const Node &node, const RectObj &aoiRect, + std::vector<std::reference_wrapper<const Feature>> &features) const { + // Does this node overlap the area of interest at all? If not, + // return without adding to the list at all. + if (!node.rect.overlaps(aoiRect)) + return; + + // Add the local features to the list. + for (const auto &pair : node.features) { + if (pair.second.overlaps(aoiRect)) { + features.push_back( + std::reference_wrapper<const Feature>(pair.first)); + } + } + + // Recurse to subnodes if they exist. + for (const auto &subnode : node.subnodes) { + search(subnode, aoiRect, features); + } + } +#endif + + void search(const Node &node, double x, double y, + std::vector<Feature> &features) const { + // Does this node overlap the area of interest at all? If not, + // return without adding to the list at all. + if (!node.rect.contains(x, y)) + return; + + // Add the local features to the list. + for (const auto &pair : node.features) { + if (pair.second.contains(x, y)) { + features.push_back(pair.first); + } + } + + // Recurse to subnodes if they exist. + for (const auto &subnode : node.subnodes) { + search(subnode, x, y, features); + } + } +}; + +} // namespace QuadTree + +NS_PROJ_END + +//! @endcond + +#endif // QUADTREE_HPP diff --git a/src/transformations/tinshift.cpp b/src/transformations/tinshift.cpp new file mode 100644 index 00000000..96e0ea4f --- /dev/null +++ b/src/transformations/tinshift.cpp @@ -0,0 +1,133 @@ +/****************************************************************************** + * Project: PROJ + * Purpose: Functionality related to triangulation transformation + * Author: Even Rouault, <even.rouault at spatialys.com> + * + ****************************************************************************** + * Copyright (c) 2020, Even Rouault, <even.rouault at spatialys.com> + * + * Permission is hereby granted, free of charge, to any person obtaining a + * copy of this software and associated documentation files (the "Software"), + * to deal in the Software without restriction, including without limitation + * the rights to use, copy, modify, merge, publish, distribute, sublicense, + * and/or sell copies of the Software, and to permit persons to whom the + * Software is furnished to do so, subject to the following conditions: + * + * The above copyright notice and this permission notice shall be included + * in all copies or substantial portions of the Software. + * + * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS + * OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, + * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL + * THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER + * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING + * FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER + * DEALINGS IN THE SOFTWARE. + *****************************************************************************/ + +#define PJ_LIB__ +#define PROJ_COMPILATION + +#include "tinshift.hpp" +#include "filemanager.hpp" +#include "proj_internal.h" + +PROJ_HEAD(tinshift, "Triangulation based transformation"); + +using namespace TINSHIFT_NAMESPACE; + +namespace { + +struct tinshiftData { + std::unique_ptr<Evaluator> evaluator{}; + + tinshiftData() = default; + + tinshiftData(const tinshiftData &) = delete; + tinshiftData &operator=(const tinshiftData &) = delete; +}; + +} // namespace + +static PJ *destructor(PJ *P, int errlev) { + if (nullptr == P) + return nullptr; + + auto Q = static_cast<struct tinshiftData *>(P->opaque); + delete Q; + P->opaque = nullptr; + + return pj_default_destructor(P, errlev); +} + +static PJ_COORD tinshift_forward_4d(PJ_COORD in, PJ *P) { + auto *Q = (struct tinshiftData *)P->opaque; + + PJ_COORD out = in; + if (!Q->evaluator->forward(in.xyz.x, in.xyz.y, in.xyz.z, out.xyz.x, + out.xyz.y, out.xyz.z)) { + return proj_coord_error(); + } + return out; +} + +static PJ_COORD tinshift_reverse_4d(PJ_COORD in, PJ *P) { + auto *Q = (struct tinshiftData *)P->opaque; + + PJ_COORD out = in; + if (!Q->evaluator->inverse(in.xyz.x, in.xyz.y, in.xyz.z, out.xyz.x, + out.xyz.y, out.xyz.z)) { + return proj_coord_error(); + } + return out; +} + +PJ *TRANSFORMATION(tinshift, 1) { + + const char *filename = pj_param(P->ctx, P->params, "sfile").s; + if (!filename) { + proj_log_error(P, "tinshift: +file= should be specified."); + return destructor(P, PJD_ERR_NO_ARGS); + } + + auto file = NS_PROJ::FileManager::open_resource_file(P->ctx, filename); + if (nullptr == file) { + proj_log_error(P, "tinshift: Cannot open %s", filename); + return destructor(P, PJD_ERR_INVALID_ARG); + } + file->seek(0, SEEK_END); + unsigned long long size = file->tell(); + // Arbitrary threshold to avoid ingesting an arbitrarily large JSON file, + // that could be a denial of service risk. 10 MB should be sufficiently + // large for any valid use ! + if (size > 10 * 1024 * 1024) { + proj_log_error(P, "tinshift: File %s too large", filename); + return destructor(P, PJD_ERR_INVALID_ARG); + } + file->seek(0); + std::string jsonStr; + jsonStr.resize(static_cast<size_t>(size)); + if (file->read(&jsonStr[0], jsonStr.size()) != jsonStr.size()) { + proj_log_error(P, "tinshift: Cannot read %s", filename); + return destructor(P, PJD_ERR_INVALID_ARG); + } + + auto Q = new tinshiftData(); + P->opaque = (void *)Q; + P->destructor = destructor; + + try { + Q->evaluator.reset(new Evaluator(TINShiftFile::parse(jsonStr))); + } catch (const std::exception &e) { + proj_log_error(P, "tinshift: invalid model: %s", e.what()); + return destructor(P, PJD_ERR_INVALID_ARG); + } + + P->destructor = destructor; + P->fwd4d = tinshift_forward_4d; + P->inv4d = tinshift_reverse_4d; + P->left = PJ_IO_UNITS_WHATEVER; + P->right = PJ_IO_UNITS_WHATEVER; + + return P; +} diff --git a/src/transformations/tinshift.hpp b/src/transformations/tinshift.hpp new file mode 100644 index 00000000..7cebab9f --- /dev/null +++ b/src/transformations/tinshift.hpp @@ -0,0 +1,257 @@ +/****************************************************************************** + * Project: PROJ + * Purpose: Functionality related to TIN based transformations + * Author: Even Rouault, <even.rouault at spatialys.com> + * + ****************************************************************************** + * Copyright (c) 2020, Even Rouault, <even.rouault at spatialys.com> + * + * Permission is hereby granted, free of charge, to any person obtaining a + * copy of this software and associated documentation files (the "Software"), + * to deal in the Software without restriction, including without limitation + * the rights to use, copy, modify, merge, publish, distribute, sublicense, + * and/or sell copies of the Software, and to permit persons to whom the + * Software is furnished to do so, subject to the following conditions: + * + * The above copyright notice and this permission notice shall be included + * in all copies or substantial portions of the Software. + * + * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS + * OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, + * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL + * THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER + * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING + * FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER + * DEALINGS IN THE SOFTWARE. + *****************************************************************************/ + +#ifndef TINSHIFT_HPP +#define TINSHIFT_HPP + +#ifdef PROJ_COMPILATION +#include "proj/internal/include_nlohmann_json.hpp" +#else +#include "nlohmann/json.hpp" +#endif + +#include <algorithm> +#include <cmath> +#include <exception> +#include <functional> +#include <limits> +#include <memory> +#include <string> +#include <vector> + +#include "quadtree.hpp" + +#ifndef TINSHIFT_NAMESPACE +#define TINSHIFT_NAMESPACE TINShift +#endif + +#include "tinshift_exceptions.hpp" + +namespace TINSHIFT_NAMESPACE { + +using json = nlohmann::json; + +// --------------------------------------------------------------------------- + +/** Content of a TINShift file. */ +class TINShiftFile { + public: + /** Parse the provided serialized JSON content and return an object. + * + * @throws ParsingException + */ + static std::unique_ptr<TINShiftFile> parse(const std::string &text); + + /** Get file type. Should always be "triangulation_file" */ + const std::string &fileType() const { return mFileType; } + + /** Get the version of the format. At time of writing, only "1.0" is known + */ + const std::string &formatVersion() const { return mFormatVersion; } + + /** Get brief descriptive name of the deformation model. */ + const std::string &name() const { return mName; } + + /** Get a string identifying the version of the deformation model. + * The format for specifying version is defined by the agency + * responsible for the deformation model. */ + const std::string &version() const { return mVersion; } + + /** Get a string identifying the license of the file. + * e.g "Create Commons Attribution 4.0 International" */ + const std::string &license() const { return mLicense; } + + /** Get a text description of the model. Intended to be longer than name() + */ + const std::string &description() const { return mDescription; } + + /** Get a text description of the model. Intended to be longer than name() + */ + const std::string &publicationDate() const { return mPublicationDate; } + + /** Basic information on the agency responsible for the model. */ + struct Authority { + std::string name{}; + std::string url{}; + std::string address{}; + std::string email{}; + }; + + /** Get basic information on the agency responsible for the model. */ + const Authority &authority() const { return mAuthority; } + + /** Hyperlink related to the model. */ + struct Link { + /** URL holding the information */ + std::string href{}; + + /** Relationship to the dataset. e.g. "about", "source", "license", + * "metadata" */ + std::string rel{}; + + /** Mime type */ + std::string type{}; + + /** Description of the link */ + std::string title{}; + }; + + /** Get links to related information. */ + const std::vector<Link> links() const { return mLinks; } + + /** Get a string identifying the CRS of source coordinates in the + * vertices. Typically "EPSG:XXXX". If the transformation is for vertical + * component, this should be the code for a compound CRS (can be + * EPSG:XXXX+YYYY where XXXX is the code of the horizontal CRS and YYYY + * the code of the vertical CRS). + * For example, for the KKJ->ETRS89 transformation, this is EPSG:2393 + * ("KKJ / Finland Uniform Coordinate System"). + * The input coordinates are assumed to + * be passed in the "normalized for visualisation" / "GIS friendly" order, + * that is longitude, latitude for geographic coordinates and + * easting, northing for projected coordinates. + * This may be empty for unspecified CRS. + */ + const std::string &inputCRS() const { return mInputCRS; } + + /** Get a string identifying the CRS of target coordinates in the + * vertices. Typically "EPSG:XXXX". If the transformation is for vertical + * component, this should be the code for a compound CRS (can be + * EPSG:XXXX+YYYY where XXXX is the code of the horizontal CRS and YYYY + * the code of the vertical CRS). + * For example, for the KKJ->ETRS89 transformation, this is EPSG:3067 + * ("ETRS89 / TM35FIN(E,N)"). + * The output coordinates will be + * returned in the "normalized for visualisation" / "GIS friendly" order, + * that is longitude, latitude for geographic coordinates and + * easting, northing for projected coordinates. + * This may be empty for unspecified CRS. + */ + const std::string &outputCRS() const { return mOutputCRS; } + + /** Return whether horizontal coordinates are transformed. */ + bool transformHorizontalComponent() const { + return mTransformHorizontalComponent; + } + + /** Return whether vertical coordinates are transformed. */ + bool transformVerticalComponent() const { + return mTransformVerticalComponent; + } + + /** Indices of vertices of a triangle */ + struct VertexIndices { + /** Index of first vertex */ + unsigned idx1; + /** Index of second vertex */ + unsigned idx2; + /** Index of third vertex */ + unsigned idx3; + }; + + /** Return number of elements per vertex of vertices() */ + unsigned verticesColumnCount() const { return mVerticesColumnCount; } + + /** Return description of triangulation vertices. + * Each vertex is described by verticesColumnCount() consecutive values. + * They are respectively: + * - the source X value + * - the source Y value + * - (if transformHorizontalComponent() is true) the target X value + * - (if transformHorizontalComponent() is true) the target Y value + * - (if transformVerticalComponent() is true) the delta Z value (to go from + * source to target Z) + * + * X is assumed to be a longitude (in degrees) or easting value. + * Y is assumed to be a latitude (in degrees) or northing value. + */ + const std::vector<double> &vertices() const { return mVertices; } + + /** Return triangles*/ + const std::vector<VertexIndices> &triangles() const { return mTriangles; } + + private: + TINShiftFile() = default; + + std::string mFileType{}; + std::string mFormatVersion{}; + std::string mName{}; + std::string mVersion{}; + std::string mLicense{}; + std::string mDescription{}; + std::string mPublicationDate{}; + Authority mAuthority{}; + std::vector<Link> mLinks{}; + std::string mInputCRS{}; + std::string mOutputCRS{}; + bool mTransformHorizontalComponent = false; + bool mTransformVerticalComponent = false; + unsigned mVerticesColumnCount = 0; + std::vector<double> mVertices{}; + std::vector<VertexIndices> mTriangles{}; +}; + +// --------------------------------------------------------------------------- + +/** Class to evaluate the transformation of a coordinate */ +class Evaluator { + public: + /** Constructor. */ + explicit Evaluator(std::unique_ptr<TINShiftFile> &&fileIn); + + /** Get file */ + const TINShiftFile &file() const { return *(mFile.get()); } + + /** Evaluate displacement of a position given by (x,y,z,t) and + * return it in (x_out,y_out_,z_out). + */ + bool forward(double x, double y, double z, double &x_out, double &y_out, + double &z_out); + + /** Apply inverse transformation. */ + bool inverse(double x, double y, double z, double &x_out, double &y_out, + double &z_out); + + private: + std::unique_ptr<TINShiftFile> mFile; + + // Reused between invokations to save memory allocations + std::vector<unsigned> mTriangleIndices{}; + + std::unique_ptr<NS_PROJ::QuadTree::QuadTree<unsigned>> mQuadTreeForward{}; + std::unique_ptr<NS_PROJ::QuadTree::QuadTree<unsigned>> mQuadTreeInverse{}; +}; + +// --------------------------------------------------------------------------- + +} // namespace TINSHIFT_NAMESPACE + +// --------------------------------------------------------------------------- + +#include "tinshift_impl.hpp" + +#endif // TINSHIFT_HPP diff --git a/src/transformations/tinshift_exceptions.hpp b/src/transformations/tinshift_exceptions.hpp new file mode 100644 index 00000000..d6d71302 --- /dev/null +++ b/src/transformations/tinshift_exceptions.hpp @@ -0,0 +1,52 @@ +/****************************************************************************** + * Project: PROJ + * Purpose: Functionality related to TIN based transformations + * Author: Even Rouault, <even.rouault at spatialys.com> + * + ****************************************************************************** + * Copyright (c) 2020, Even Rouault, <even.rouault at spatialys.com> + * + * Permission is hereby granted, free of charge, to any person obtaining a + * copy of this software and associated documentation files (the "Software"), + * to deal in the Software without restriction, including without limitation + * the rights to use, copy, modify, merge, publish, distribute, sublicense, + * and/or sell copies of the Software, and to permit persons to whom the + * Software is furnished to do so, subject to the following conditions: + * + * The above copyright notice and this permission notice shall be included + * in all copies or substantial portions of the Software. + * + * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS + * OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, + * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL + * THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER + * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING + * FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER + * DEALINGS IN THE SOFTWARE. + *****************************************************************************/ + +#ifndef TINSHIFT_NAMESPACE +#error "Should be included only by tinshift.hpp" +#endif + +#include <exception> + +namespace TINSHIFT_NAMESPACE { + +// --------------------------------------------------------------------------- + +/** Parsing exception. */ +class ParsingException : public std::exception { + public: + explicit ParsingException(const std::string &msg) : msg_(msg) {} + const char *what() const noexcept override; + + private: + std::string msg_; +}; + +const char *ParsingException::what() const noexcept { return msg_.c_str(); } + +// --------------------------------------------------------------------------- + +} // namespace DEFORMATON_MODEL_NAMESPACE diff --git a/src/transformations/tinshift_impl.hpp b/src/transformations/tinshift_impl.hpp new file mode 100644 index 00000000..90f2b2f8 --- /dev/null +++ b/src/transformations/tinshift_impl.hpp @@ -0,0 +1,555 @@ +/****************************************************************************** + * Project: PROJ + * Purpose: Functionality related to TIN based transformations + * Author: Even Rouault, <even.rouault at spatialys.com> + * + ****************************************************************************** + * Copyright (c) 2020, Even Rouault, <even.rouault at spatialys.com> + * + * Permission is hereby granted, free of charge, to any person obtaining a + * copy of this software and associated documentation files (the "Software"), + * to deal in the Software without restriction, including without limitation + * the rights to use, copy, modify, merge, publish, distribute, sublicense, + * and/or sell copies of the Software, and to permit persons to whom the + * Software is furnished to do so, subject to the following conditions: + * + * The above copyright notice and this permission notice shall be included + * in all copies or substantial portions of the Software. + * + * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS + * OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, + * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL + * THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER + * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING + * FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER + * DEALINGS IN THE SOFTWARE. + *****************************************************************************/ + +#ifndef TINSHIFT_NAMESPACE +#error "Should be included only by tinshift.hpp" +#endif + +#include <limits> + +namespace TINSHIFT_NAMESPACE { + +// --------------------------------------------------------------------------- + +static std::string getString(const json &j, const char *key, bool optional) { + if (!j.contains(key)) { + if (optional) { + return std::string(); + } + throw ParsingException(std::string("Missing \"") + key + "\" key"); + } + const json v = j[key]; + if (!v.is_string()) { + throw ParsingException(std::string("The value of \"") + key + + "\" should be a string"); + } + return v.get<std::string>(); +} + +static std::string getReqString(const json &j, const char *key) { + return getString(j, key, false); +} + +static std::string getOptString(const json &j, const char *key) { + return getString(j, key, true); +} + +// --------------------------------------------------------------------------- + +static json getArrayMember(const json &j, const char *key) { + if (!j.contains(key)) { + throw ParsingException(std::string("Missing \"") + key + "\" key"); + } + const json obj = j[key]; + if (!obj.is_array()) { + throw ParsingException(std::string("The value of \"") + key + + "\" should be a array"); + } + return obj; +} + +// --------------------------------------------------------------------------- + +std::unique_ptr<TINShiftFile> TINShiftFile::parse(const std::string &text) { + std::unique_ptr<TINShiftFile> tinshiftFile(new TINShiftFile()); + json j; + try { + j = json::parse(text); + } catch (const std::exception &e) { + throw ParsingException(e.what()); + } + if (!j.is_object()) { + throw ParsingException("Not an object"); + } + tinshiftFile->mFileType = getReqString(j, "file_type"); + tinshiftFile->mFormatVersion = getReqString(j, "format_version"); + tinshiftFile->mName = getOptString(j, "name"); + tinshiftFile->mVersion = getOptString(j, "version"); + tinshiftFile->mLicense = getOptString(j, "license"); + tinshiftFile->mDescription = getOptString(j, "description"); + tinshiftFile->mPublicationDate = getOptString(j, "publication_date"); + + if (j.contains("authority")) { + const json jAuthority = j["authority"]; + if (!jAuthority.is_object()) { + throw ParsingException("authority is not a object"); + } + tinshiftFile->mAuthority.name = getOptString(jAuthority, "name"); + tinshiftFile->mAuthority.url = getOptString(jAuthority, "url"); + tinshiftFile->mAuthority.address = getOptString(jAuthority, "address"); + tinshiftFile->mAuthority.email = getOptString(jAuthority, "email"); + } + + if (j.contains("links")) { + const json jLinks = j["links"]; + if (!jLinks.is_array()) { + throw ParsingException("links is not an array"); + } + for (const json &jLink : jLinks) { + if (!jLink.is_object()) { + throw ParsingException("links[] item is not an object"); + } + Link link; + link.href = getOptString(jLink, "href"); + link.rel = getOptString(jLink, "rel"); + link.type = getOptString(jLink, "type"); + link.title = getOptString(jLink, "title"); + tinshiftFile->mLinks.emplace_back(std::move(link)); + } + } + tinshiftFile->mInputCRS = getOptString(j, "input_crs"); + tinshiftFile->mOutputCRS = getOptString(j, "output_crs"); + + const auto jTransformedComponents = + getArrayMember(j, "transformed_components"); + for (const json &jComp : jTransformedComponents) { + if (!jComp.is_string()) { + throw ParsingException( + "transformed_components[] item is not a string"); + } + const auto jCompStr = jComp.get<std::string>(); + if (jCompStr == "horizontal") { + tinshiftFile->mTransformHorizontalComponent = true; + } else if (jCompStr == "vertical") { + tinshiftFile->mTransformVerticalComponent = true; + } else { + throw ParsingException("transformed_components[] = " + jCompStr + + " is not handled"); + } + } + + const auto jVerticesColumns = getArrayMember(j, "vertices_columns"); + int sourceXCol = -1; + int sourceYCol = -1; + int sourceZCol = -1; + int targetXCol = -1; + int targetYCol = -1; + int targetZCol = -1; + int offsetZCol = -1; + for (size_t i = 0; i < jVerticesColumns.size(); ++i) { + const json &jColumn = jVerticesColumns[i]; + if (!jColumn.is_string()) { + throw ParsingException("vertices_columns[] item is not a string"); + } + const auto jColumnStr = jColumn.get<std::string>(); + if (jColumnStr == "source_x") { + sourceXCol = static_cast<int>(i); + } else if (jColumnStr == "source_y") { + sourceYCol = static_cast<int>(i); + } else if (jColumnStr == "source_z") { + sourceZCol = static_cast<int>(i); + } else if (jColumnStr == "target_x") { + targetXCol = static_cast<int>(i); + } else if (jColumnStr == "target_y") { + targetYCol = static_cast<int>(i); + } else if (jColumnStr == "target_z") { + targetZCol = static_cast<int>(i); + } else if (jColumnStr == "offset_z") { + offsetZCol = static_cast<int>(i); + } + } + if (sourceXCol < 0) { + throw ParsingException( + "source_x must be specified in vertices_columns[]"); + } + if (sourceYCol < 0) { + throw ParsingException( + "source_y must be specified in vertices_columns[]"); + } + if (tinshiftFile->mTransformHorizontalComponent) { + if (targetXCol < 0) { + throw ParsingException( + "target_x must be specified in vertices_columns[]"); + } + if (targetYCol < 0) { + throw ParsingException( + "target_y must be specified in vertices_columns[]"); + } + } + if (tinshiftFile->mTransformVerticalComponent) { + if (offsetZCol >= 0) { + // do nothing + } else { + if (sourceZCol < 0) { + throw ParsingException("source_z or delta_z must be specified " + "in vertices_columns[]"); + } + if (targetZCol < 0) { + throw ParsingException( + "target_z must be specified in vertices_columns[]"); + } + } + } + + const auto jTrianglesColumns = getArrayMember(j, "triangles_columns"); + int idxVertex1Col = -1; + int idxVertex2Col = -1; + int idxVertex3Col = -1; + for (size_t i = 0; i < jTrianglesColumns.size(); ++i) { + const json &jColumn = jTrianglesColumns[i]; + if (!jColumn.is_string()) { + throw ParsingException("triangles_columns[] item is not a string"); + } + const auto jColumnStr = jColumn.get<std::string>(); + if (jColumnStr == "idx_vertex1") { + idxVertex1Col = static_cast<int>(i); + } else if (jColumnStr == "idx_vertex2") { + idxVertex2Col = static_cast<int>(i); + } else if (jColumnStr == "idx_vertex3") { + idxVertex3Col = static_cast<int>(i); + } + } + if (idxVertex1Col < 0) { + throw ParsingException( + "idx_vertex1 must be specified in triangles_columns[]"); + } + if (idxVertex2Col < 0) { + throw ParsingException( + "idx_vertex2 must be specified in triangles_columns[]"); + } + if (idxVertex3Col < 0) { + throw ParsingException( + "idx_vertex3 must be specified in triangles_columns[]"); + } + + const auto jVertices = getArrayMember(j, "vertices"); + tinshiftFile->mVerticesColumnCount = 2; + if (tinshiftFile->mTransformHorizontalComponent) + tinshiftFile->mVerticesColumnCount += 2; + if (tinshiftFile->mTransformVerticalComponent) + tinshiftFile->mVerticesColumnCount += 1; + + tinshiftFile->mVertices.reserve(tinshiftFile->mVerticesColumnCount * + jVertices.size()); + for (const auto &jVertex : jVertices) { + if (!jVertex.is_array()) { + throw ParsingException("vertices[] item is not an array"); + } + if (jVertex.size() != jVerticesColumns.size()) { + throw ParsingException( + "vertices[] item has not expected number of elements"); + } + if (!jVertex[sourceXCol].is_number()) { + throw ParsingException("vertices[][] item is not a number"); + } + tinshiftFile->mVertices.push_back(jVertex[sourceXCol].get<double>()); + if (!jVertex[sourceYCol].is_number()) { + throw ParsingException("vertices[][] item is not a number"); + } + tinshiftFile->mVertices.push_back(jVertex[sourceYCol].get<double>()); + if (tinshiftFile->mTransformHorizontalComponent) { + if (!jVertex[targetXCol].is_number()) { + throw ParsingException("vertices[][] item is not a number"); + } + tinshiftFile->mVertices.push_back( + jVertex[targetXCol].get<double>()); + if (!jVertex[targetYCol].is_number()) { + throw ParsingException("vertices[][] item is not a number"); + } + tinshiftFile->mVertices.push_back( + jVertex[targetYCol].get<double>()); + } + if (tinshiftFile->mTransformVerticalComponent) { + if (offsetZCol >= 0) { + if (!jVertex[offsetZCol].is_number()) { + throw ParsingException("vertices[][] item is not a number"); + } + tinshiftFile->mVertices.push_back( + jVertex[offsetZCol].get<double>()); + } else { + if (!jVertex[sourceZCol].is_number()) { + throw ParsingException("vertices[][] item is not a number"); + } + const double sourceZ = jVertex[sourceZCol].get<double>(); + if (!jVertex[targetZCol].is_number()) { + throw ParsingException("vertices[][] item is not a number"); + } + const double targetZ = jVertex[targetZCol].get<double>(); + tinshiftFile->mVertices.push_back(targetZ - sourceZ); + } + } + } + + const auto jTriangles = getArrayMember(j, "triangles"); + tinshiftFile->mTriangles.reserve(jTriangles.size()); + for (const auto &jTriangle : jTriangles) { + if (!jTriangle.is_array()) { + throw ParsingException("triangles[] item is not an array"); + } + if (jTriangle.size() != jTrianglesColumns.size()) { + throw ParsingException( + "triangles[] item has not expected number of elements"); + } + + if (jTriangle[idxVertex1Col].type() != json::value_t::number_unsigned) { + throw ParsingException("triangles[][] item is not an integer"); + } + const unsigned vertex1 = jTriangle[idxVertex1Col].get<unsigned>(); + if (vertex1 >= jVertices.size()) { + throw ParsingException("Invalid value for a vertex index"); + } + + if (jTriangle[idxVertex2Col].type() != json::value_t::number_unsigned) { + throw ParsingException("triangles[][] item is not an integer"); + } + const unsigned vertex2 = jTriangle[idxVertex2Col].get<unsigned>(); + if (vertex2 >= jVertices.size()) { + throw ParsingException("Invalid value for a vertex index"); + } + + if (jTriangle[idxVertex3Col].type() != json::value_t::number_unsigned) { + throw ParsingException("triangles[][] item is not an integer"); + } + const unsigned vertex3 = jTriangle[idxVertex3Col].get<unsigned>(); + if (vertex3 >= jVertices.size()) { + throw ParsingException("Invalid value for a vertex index"); + } + + VertexIndices vi; + vi.idx1 = vertex1; + vi.idx2 = vertex2; + vi.idx3 = vertex3; + tinshiftFile->mTriangles.push_back(vi); + } + + return tinshiftFile; +} + +// --------------------------------------------------------------------------- + +static NS_PROJ::QuadTree::RectObj GetBounds(const TINShiftFile &file, + bool forward) { + NS_PROJ::QuadTree::RectObj rect; + rect.minx = std::numeric_limits<double>::max(); + rect.miny = std::numeric_limits<double>::max(); + rect.maxx = -std::numeric_limits<double>::max(); + rect.maxy = -std::numeric_limits<double>::max(); + const auto &vertices = file.vertices(); + const unsigned colCount = file.verticesColumnCount(); + const int idxX = file.transformHorizontalComponent() && !forward ? 2 : 0; + const int idxY = file.transformHorizontalComponent() && !forward ? 3 : 1; + for (size_t i = 0; i + colCount - 1 < vertices.size(); i += colCount) { + const double x = vertices[i + idxX]; + const double y = vertices[i + idxY]; + rect.minx = std::min(rect.minx, x); + rect.miny = std::min(rect.miny, y); + rect.maxx = std::max(rect.maxx, x); + rect.maxy = std::max(rect.maxy, y); + } + return rect; +} + +// --------------------------------------------------------------------------- + +static std::unique_ptr<NS_PROJ::QuadTree::QuadTree<unsigned>> +BuildQuadTree(const TINShiftFile &file, bool forward) { + auto quadtree = std::unique_ptr<NS_PROJ::QuadTree::QuadTree<unsigned>>( + new NS_PROJ::QuadTree::QuadTree<unsigned>(GetBounds(file, forward))); + const auto &triangles = file.triangles(); + const auto &vertices = file.vertices(); + const int idxX = file.transformHorizontalComponent() && !forward ? 2 : 0; + const int idxY = file.transformHorizontalComponent() && !forward ? 3 : 1; + const unsigned colCount = file.verticesColumnCount(); + for (size_t i = 0; i < triangles.size(); ++i) { + const unsigned i1 = triangles[i].idx1; + const unsigned i2 = triangles[i].idx2; + const unsigned i3 = triangles[i].idx3; + const double x1 = vertices[i1 * colCount + idxX]; + const double y1 = vertices[i1 * colCount + idxY]; + const double x2 = vertices[i2 * colCount + idxX]; + const double y2 = vertices[i2 * colCount + idxY]; + const double x3 = vertices[i3 * colCount + idxX]; + const double y3 = vertices[i3 * colCount + idxY]; + NS_PROJ::QuadTree::RectObj rect; + rect.minx = x1; + rect.miny = y1; + rect.maxx = x1; + rect.maxy = y1; + rect.minx = std::min(rect.minx, x2); + rect.miny = std::min(rect.miny, y2); + rect.maxx = std::max(rect.maxx, x2); + rect.maxy = std::max(rect.maxy, y2); + rect.minx = std::min(rect.minx, x3); + rect.miny = std::min(rect.miny, y3); + rect.maxx = std::max(rect.maxx, x3); + rect.maxy = std::max(rect.maxy, y3); + quadtree->insert(static_cast<unsigned>(i), rect); + } + + return quadtree; +} + +// --------------------------------------------------------------------------- + +Evaluator::Evaluator(std::unique_ptr<TINShiftFile> &&fileIn) + : mFile(std::move(fileIn)) {} + +// --------------------------------------------------------------------------- + +static const TINShiftFile::VertexIndices * +FindTriangle(const TINShiftFile &file, + const NS_PROJ::QuadTree::QuadTree<unsigned> &quadtree, + std::vector<unsigned> &triangleIndices, double x, double y, + bool forward, double &lambda1, double &lambda2, double &lambda3) { +#define USE_QUADTREE +#ifdef USE_QUADTREE + triangleIndices.clear(); + quadtree.search(x, y, triangleIndices); +#endif + const auto &triangles = file.triangles(); + const auto &vertices = file.vertices(); + constexpr double EPS = 1e-10; + const int idxX = file.transformHorizontalComponent() && !forward ? 2 : 0; + const int idxY = file.transformHorizontalComponent() && !forward ? 3 : 1; + const unsigned colCount = file.verticesColumnCount(); +#ifdef USE_QUADTREE + for (unsigned i : triangleIndices) +#else + for (size_t i = 0; i < triangles.size(); ++i) +#endif + { + const auto &triangle = triangles[i]; + const unsigned i1 = triangle.idx1; + const unsigned i2 = triangle.idx2; + const unsigned i3 = triangle.idx3; + const double x1 = vertices[i1 * colCount + idxX]; + const double y1 = vertices[i1 * colCount + idxY]; + const double x2 = vertices[i2 * colCount + idxX]; + const double y2 = vertices[i2 * colCount + idxY]; + const double x3 = vertices[i3 * colCount + idxX]; + const double y3 = vertices[i3 * colCount + idxY]; + const double det_T = (y2 - y3) * (x1 - x3) + (x3 - x2) * (y1 - y3); + lambda1 = ((y2 - y3) * (x - x3) + (x3 - x2) * (y - y3)) / det_T; + lambda2 = ((y3 - y1) * (x - x3) + (x1 - x3) * (y - y3)) / det_T; + if (lambda1 >= -EPS && lambda1 <= 1 + EPS && lambda2 >= -EPS && + lambda2 <= 1 + EPS) { + lambda3 = 1 - lambda1 - lambda2; + if (lambda3 >= 0) { + return ▵ + } + } + } + return nullptr; +} + +// --------------------------------------------------------------------------- + +bool Evaluator::forward(double x, double y, double z, double &x_out, + double &y_out, double &z_out) { + if (!mQuadTreeForward) + mQuadTreeForward = BuildQuadTree(*(mFile.get()), true); + + double lambda1 = 0.0; + double lambda2 = 0.0; + double lambda3 = 0.0; + const auto *triangle = + FindTriangle(*mFile, *mQuadTreeForward, mTriangleIndices, x, y, true, + lambda1, lambda2, lambda3); + if (!triangle) + return false; + const auto &vertices = mFile->vertices(); + const unsigned i1 = triangle->idx1; + const unsigned i2 = triangle->idx2; + const unsigned i3 = triangle->idx3; + const unsigned colCount = mFile->verticesColumnCount(); + if (mFile->transformHorizontalComponent()) { + constexpr unsigned idxTargetColX = 2; + constexpr unsigned idxTargetColY = 3; + x_out = vertices[i1 * colCount + idxTargetColX] * lambda1 + + vertices[i2 * colCount + idxTargetColX] * lambda2 + + vertices[i3 * colCount + idxTargetColX] * lambda3; + y_out = vertices[i1 * colCount + idxTargetColY] * lambda1 + + vertices[i2 * colCount + idxTargetColY] * lambda2 + + vertices[i3 * colCount + idxTargetColY] * lambda3; + } else { + x_out = x; + y_out = y; + } + if (mFile->transformVerticalComponent()) { + const int idxCol = mFile->transformHorizontalComponent() ? 4 : 2; + z_out = z + (vertices[i1 * colCount + idxCol] * lambda1 + + vertices[i2 * colCount + idxCol] * lambda2 + + vertices[i3 * colCount + idxCol] * lambda3); + } else { + z_out = z; + } + return true; +} + +// --------------------------------------------------------------------------- + +bool Evaluator::inverse(double x, double y, double z, double &x_out, + double &y_out, double &z_out) { + NS_PROJ::QuadTree::QuadTree<unsigned> *quadtree; + if (!mFile->transformHorizontalComponent() && + mFile->transformVerticalComponent()) { + if (!mQuadTreeForward) + mQuadTreeForward = BuildQuadTree(*(mFile.get()), true); + quadtree = mQuadTreeForward.get(); + } else { + if (!mQuadTreeInverse) + mQuadTreeInverse = BuildQuadTree(*(mFile.get()), false); + quadtree = mQuadTreeInverse.get(); + } + + double lambda1 = 0.0; + double lambda2 = 0.0; + double lambda3 = 0.0; + const auto *triangle = FindTriangle(*mFile, *quadtree, mTriangleIndices, x, + y, false, lambda1, lambda2, lambda3); + if (!triangle) + return false; + const auto &vertices = mFile->vertices(); + const unsigned i1 = triangle->idx1; + const unsigned i2 = triangle->idx2; + const unsigned i3 = triangle->idx3; + const unsigned colCount = mFile->verticesColumnCount(); + if (mFile->transformHorizontalComponent()) { + constexpr unsigned idxTargetColX = 0; + constexpr unsigned idxTargetColY = 1; + x_out = vertices[i1 * colCount + idxTargetColX] * lambda1 + + vertices[i2 * colCount + idxTargetColX] * lambda2 + + vertices[i3 * colCount + idxTargetColX] * lambda3; + y_out = vertices[i1 * colCount + idxTargetColY] * lambda1 + + vertices[i2 * colCount + idxTargetColY] * lambda2 + + vertices[i3 * colCount + idxTargetColY] * lambda3; + } else { + x_out = x; + y_out = y; + } + if (mFile->transformVerticalComponent()) { + const int idxCol = mFile->transformHorizontalComponent() ? 4 : 2; + z_out = z - (vertices[i1 * colCount + idxCol] * lambda1 + + vertices[i2 * colCount + idxCol] * lambda2 + + vertices[i3 * colCount + idxCol] * lambda3); + } else { + z_out = z; + } + return true; +} + +} // namespace TINSHIFT_NAMESPACE diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt index 9fc29103..c60f344b 100644 --- a/test/CMakeLists.txt +++ b/test/CMakeLists.txt @@ -17,6 +17,7 @@ proj_add_gie_test("adams_ws2" "gie/adams_ws2.gie") proj_add_gie_test("guyou" "gie/guyou.gie") proj_add_gie_test("peirce_q" "gie/peirce_q.gie") proj_add_gie_test("defmodel" "gie/defmodel.gie") +proj_add_gie_test("tinshift" "gie/tinshift.gie") # GIGS tests. Uncommented tests are expected to fail due to issues with # various projections. Should be investigated further and fixed. diff --git a/test/gie/Makefile.am b/test/gie/Makefile.am index e6bb9582..20f0ec8f 100644 --- a/test/gie/Makefile.am +++ b/test/gie/Makefile.am @@ -16,7 +16,8 @@ EXTRA_DIST = 4D-API_cs2cs-style.gie \ adams_ws2.gie \ guyou.gie \ peirce_q.gie \ - defmodel.gie + defmodel.gie \ + tinshift.gie PROJ_LIB ?= ../../data/for_tests @@ -68,4 +69,7 @@ peirce_q: peirce_q.gie defmodel: defmodel.gie PROJ_SKIP_READ_USER_WRITABLE_DIRECTORY=YES PROJ_LIB=$(PROJ_LIB) $(GIEEXE) $< -check-local: 4D-API-cs2cs-style GDA axisswap builtins deformation ellipsoid more_builtins unitconvert DHDN_ETRS89 geotiff_grids adams_hemi adams_ws1 adams_ws2 guyou peirce_q defmodel +tinshift: tinshift.gie + PROJ_SKIP_READ_USER_WRITABLE_DIRECTORY=YES PROJ_LIB=$(PROJ_LIB) $(GIEEXE) $< + +check-local: 4D-API-cs2cs-style GDA axisswap builtins deformation ellipsoid more_builtins unitconvert DHDN_ETRS89 geotiff_grids adams_hemi adams_ws1 adams_ws2 guyou peirce_q defmodel tinshift diff --git a/test/gie/tinshift.gie b/test/gie/tinshift.gie new file mode 100644 index 00000000..dff62693 --- /dev/null +++ b/test/gie/tinshift.gie @@ -0,0 +1,50 @@ + +------------------------------------------------------------------------------- +=============================================================================== +Test +proj=tinshift +=============================================================================== + +<gie-strict> + +# Missing +file +operation +proj=tinshift +expect failure errno no_args + +# +file doesn't point to an existing file +operation +proj=tinshift +file=i_do_not_exist +expect failure errno invalid_arg + +# Not a JSON file +operation +proj=tinshift +file=proj.ini +expect failure errno invalid_arg + + +# Tests on a file without explicit CRS +operation +proj=tinshift +file=tests/tinshift_crs_implicit.json +accept 2 49 +expect 2.1 49.1 +roundtrip 1 + +accept 0 0 +expect failure + +direction inverse +accept 0 0 +expect failure + + +# Tests on a file with explicit CRS +operation +proj=tinshift +file=tests/tinshift_simplified_kkj_etrs.json +tolerance 0.1 mm +# Verified with https://kartta.paikkatietoikkuna.fi/?lang=en with EPSG:2393 to EPSG:3067 +accept 3210000.0000 6700000.0000 +expect 209948.3217 6697187.0009 +roundtrip 1 + +operation +proj=tinshift +file=tests/tinshift_simplified_n60_n2000.json +tolerance 0.1 mm +accept 3210000.0000 6700000.0000 10.0 +expect 3210000.0000 6700000.0000 10.2886 +roundtrip 1 + +</gie-strict> diff --git a/test/unit/CMakeLists.txt b/test/unit/CMakeLists.txt index 04c6acb6..2c0c19a9 100644 --- a/test/unit/CMakeLists.txt +++ b/test/unit/CMakeLists.txt @@ -181,3 +181,14 @@ target_link_libraries(test_defmodel add_test(NAME test_defmodel COMMAND test_defmodel) set_property(TEST test_defmodel PROPERTY ENVIRONMENT ${PROJ_TEST_ENVIRONMENT}) + + +add_executable(test_tinshift + main.cpp + test_tinshift.cpp) +target_link_libraries(test_tinshift + GTest::gtest + ${PROJ_LIBRARIES}) +add_test(NAME test_tinshift COMMAND test_tinshift) +set_property(TEST test_tinshift + PROPERTY ENVIRONMENT ${PROJ_TEST_ENVIRONMENT}) diff --git a/test/unit/Makefile.am b/test/unit/Makefile.am index 2dd0dbff..483cb0bd 100644 --- a/test/unit/Makefile.am +++ b/test/unit/Makefile.am @@ -19,6 +19,7 @@ noinst_PROGRAMS += gie_self_tests noinst_PROGRAMS += include_proj_h_from_c noinst_PROGRAMS += test_network noinst_PROGRAMS += test_defmodel +noinst_PROGRAMS += test_tinshift pj_transform_test_SOURCES = pj_transform_test.cpp main.cpp pj_transform_test_LDADD = ../../src/libproj.la @GTEST_LIBS@ @@ -77,6 +78,12 @@ test_defmodel_LDADD = ../../src/libproj.la @GTEST_LIBS@ test_defmodel-check: test_defmodel PROJ_LIB=$(PROJ_LIB) ./test_defmodel +test_tinshift_SOURCES = test_tinshift.cpp main.cpp +test_tinshift_LDADD = ../../src/libproj.la @GTEST_LIBS@ + +test_tinshift-check: test_tinshift + PROJ_LIB=$(PROJ_LIB) ./test_tinshift + check-local: pj_transform_test-check pj_phi2_test-check proj_errno_string_test-check \ proj_angular_io_test-check proj_context_test-check test_cpp_api-check \ - gie_self_tests-check test_network-check test_defmodel-check + gie_self_tests-check test_network-check test_defmodel-check test_tinshift-check diff --git a/test/unit/test_tinshift.cpp b/test/unit/test_tinshift.cpp new file mode 100644 index 00000000..96da6e25 --- /dev/null +++ b/test/unit/test_tinshift.cpp @@ -0,0 +1,211 @@ +/****************************************************************************** + * + * Project: PROJ + * Purpose: Test TIN shift + * Author: Even Rouault <even dot rouault at spatialys dot com> + * + ****************************************************************************** + * Copyright (c) 2020, Even Rouault <even dot rouault at spatialys dot com> + * + * Permission is hereby granted, free of charge, to any person obtaining a + * copy of this software and associated documentation files (the "Software"), + * to deal in the Software without restriction, including without limitation + * the rights to use, copy, modify, merge, publish, distribute, sublicense, + * and/or sell copies of the Software, and to permit persons to whom the + * Software is furnished to do so, subject to the following conditions: + * + * The above copyright notice and this permission notice shall be included + * in all copies or substantial portions of the Software. + * + * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS + * OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, + * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL + * THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER + * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING + * FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER + * DEALINGS IN THE SOFTWARE. + ****************************************************************************/ + +#include "gtest_include.h" + +#define PROJ_COMPILATION +#define TINSHIFT_NAMESPACE TestTINShift +#include "transformations/tinshift.hpp" + +using namespace TINSHIFT_NAMESPACE; + +namespace { + +static json getMinValidContent() { + json j; + j["file_type"] = "triangulation_file"; + j["format_version"] = "1.0"; + j["input_crs"] = "EPSG:2393"; + j["output_crs"] = "EPSG:3067"; + j["transformed_components"] = {"horizontal"}; + j["vertices_columns"] = {"source_x", "source_y", "target_x", "target_y"}; + j["triangles_columns"] = {"idx_vertex1", "idx_vertex2", "idx_vertex3"}; + j["vertices"] = {{0, 0, 101, 101}, {0, 1, 100, 101}, {1, 1, 100, 100}}; + j["triangles"] = {{0, 1, 2}}; + + return j; +} + +// --------------------------------------------------------------------------- + +TEST(tinshift, basic) { + EXPECT_THROW(TINShiftFile::parse("foo"), ParsingException); + EXPECT_THROW(TINShiftFile::parse("null"), ParsingException); + EXPECT_THROW(TINShiftFile::parse("{}"), ParsingException); + + const auto jMinValid(getMinValidContent()); + { + auto f = TINShiftFile::parse(jMinValid.dump()); + EXPECT_EQ(f->fileType(), "triangulation_file"); + EXPECT_EQ(f->formatVersion(), "1.0"); + EXPECT_EQ(f->inputCRS(), "EPSG:2393"); + EXPECT_EQ(f->outputCRS(), "EPSG:3067"); + + auto eval = Evaluator(std::move(f)); + double x_out = 0; + double y_out = 0; + double z_out = 0; + + EXPECT_FALSE(eval.forward(-0.1, 0.0, 1000.0, x_out, y_out, z_out)); + + EXPECT_TRUE(eval.forward(0.0, 0.0, 1000.0, x_out, y_out, z_out)); + EXPECT_EQ(x_out, 101.0); + EXPECT_EQ(y_out, 101.0); + EXPECT_EQ(z_out, 1000.0); + + EXPECT_TRUE(eval.forward(0.0, 1.0, 1000.0, x_out, y_out, z_out)); + EXPECT_EQ(x_out, 100.0); + EXPECT_EQ(y_out, 101.0); + EXPECT_EQ(z_out, 1000.0); + + EXPECT_TRUE(eval.forward(1.0, 1.0, 1000.0, x_out, y_out, z_out)); + EXPECT_EQ(x_out, 100.0); + EXPECT_EQ(y_out, 100.0); + EXPECT_EQ(z_out, 1000.0); + + EXPECT_TRUE(eval.forward(0.0, 0.5, 1000.0, x_out, y_out, z_out)); + EXPECT_EQ(x_out, 100.5); + EXPECT_EQ(y_out, 101.0); + EXPECT_EQ(z_out, 1000.0); + + EXPECT_TRUE(eval.forward(0.5, 0.5, 1000.0, x_out, y_out, z_out)); + EXPECT_EQ(x_out, 100.5); + EXPECT_EQ(y_out, 100.5); + EXPECT_EQ(z_out, 1000.0); + + EXPECT_TRUE(eval.forward(0.5, 0.75, 1000.0, x_out, y_out, z_out)); + EXPECT_EQ(x_out, 100.25); + EXPECT_EQ(y_out, 100.5); + EXPECT_EQ(z_out, 1000.0); + + EXPECT_TRUE(eval.inverse(100.25, 100.5, 1000.0, x_out, y_out, z_out)); + EXPECT_EQ(x_out, 0.5); + EXPECT_EQ(y_out, 0.75); + EXPECT_EQ(z_out, 1000.0); + } + + // Vertical only, through source_z / target_z + { + auto j(jMinValid); + j["transformed_components"] = {"vertical"}; + j["vertices_columns"] = {"source_x", "source_y", "source_z", + "target_z"}; + j["triangles_columns"] = {"idx_vertex1", "idx_vertex2", "idx_vertex3"}; + j["vertices"] = { + {0, 0, 10.5, 10.6}, {0, 1, 15.0, 15.2}, {1, 1, 17.5, 18.0}}; + j["triangles"] = {{0, 1, 2}}; + + auto f = TINShiftFile::parse(j.dump()); + auto eval = Evaluator(std::move(f)); + double x_out = 0; + double y_out = 0; + double z_out = 0; + + EXPECT_TRUE(eval.forward(0.0, 0.0, 1000.0, x_out, y_out, z_out)); + EXPECT_EQ(x_out, 0.0); + EXPECT_EQ(y_out, 0.0); + EXPECT_EQ(z_out, 1000.1); + + EXPECT_TRUE(eval.forward(0.5, 0.75, 1000.0, x_out, y_out, z_out)); + EXPECT_EQ(x_out, 0.5); + EXPECT_EQ(y_out, 0.75); + EXPECT_EQ(z_out, 1000.325); + + EXPECT_TRUE(eval.inverse(0.5, 0.75, 1000.325, x_out, y_out, z_out)); + EXPECT_EQ(x_out, 0.5); + EXPECT_EQ(y_out, 0.75); + EXPECT_EQ(z_out, 1000.0); + } + + // Vertical only, through offset_z + { + auto j(jMinValid); + j["transformed_components"] = {"vertical"}; + j["vertices_columns"] = {"source_x", "source_y", "offset_z"}; + j["triangles_columns"] = {"idx_vertex1", "idx_vertex2", "idx_vertex3"}; + j["vertices"] = {{0, 0, 0.1}, {0, 1, 0.2}, {1, 1, 0.5}}; + j["triangles"] = {{0, 1, 2}}; + + auto f = TINShiftFile::parse(j.dump()); + auto eval = Evaluator(std::move(f)); + double x_out = 0; + double y_out = 0; + double z_out = 0; + + EXPECT_TRUE(eval.forward(0.0, 0.0, 1000.0, x_out, y_out, z_out)); + EXPECT_EQ(x_out, 0.0); + EXPECT_EQ(y_out, 0.0); + EXPECT_EQ(z_out, 1000.1); + + EXPECT_TRUE(eval.forward(0.5, 0.75, 1000.0, x_out, y_out, z_out)); + EXPECT_EQ(x_out, 0.5); + EXPECT_EQ(y_out, 0.75); + EXPECT_EQ(z_out, 1000.325); + + EXPECT_TRUE(eval.inverse(0.5, 0.75, 1000.325, x_out, y_out, z_out)); + EXPECT_EQ(x_out, 0.5); + EXPECT_EQ(y_out, 0.75); + EXPECT_EQ(z_out, 1000.0); + } + + // Horizontal + vertical + { + auto j(jMinValid); + j["transformed_components"] = {"horizontal", "vertical"}; + j["vertices_columns"] = {"source_x", "source_y", "target_x", "target_y", + "offset_z"}; + j["triangles_columns"] = {"idx_vertex1", "idx_vertex2", "idx_vertex3"}; + j["vertices"] = {{0, 0, 101, 101, 0.1}, + {0, 1, 100, 101, 0.2}, + {1, 1, 100, 100, 0.5}}; + j["triangles"] = {{0, 1, 2}}; + + auto f = TINShiftFile::parse(j.dump()); + auto eval = Evaluator(std::move(f)); + double x_out = 0; + double y_out = 0; + double z_out = 0; + + EXPECT_TRUE(eval.forward(0.0, 0.0, 1000.0, x_out, y_out, z_out)); + EXPECT_EQ(x_out, 101.0); + EXPECT_EQ(y_out, 101.0); + EXPECT_EQ(z_out, 1000.1); + + EXPECT_TRUE(eval.forward(0.5, 0.75, 1000.0, x_out, y_out, z_out)); + EXPECT_EQ(x_out, 100.25); + EXPECT_EQ(y_out, 100.5); + EXPECT_EQ(z_out, 1000.325); + + EXPECT_TRUE(eval.inverse(100.25, 100.5, 1000.325, x_out, y_out, z_out)); + EXPECT_EQ(x_out, 0.5); + EXPECT_EQ(y_out, 0.75); + EXPECT_EQ(z_out, 1000.0); + } +} + +} // namespace |