Subscribe to this thread
Home - Cutting Edge / All posts - Manifold System

7,972 post(s)
#19-Jun-18 15:29

SHA256: cffa474acbbea37515c15b0359337177f0bfe4beb85f917780e054bc84651d75

SHA256: e35b5dc2ed82cf8cfba2e1c67f2dadc43e724cf33d66f205a485361e7776d932


7,972 post(s)
#19-Jun-18 15:31

Interpolating Surfaces

There are query functions for interpolating surfaces.

Two caveats: First, the performance is not final. We have multiple big improvements currently in testing, they make for a big difference. Second, the Transform pane does not yet include point-and-click templates for interpolating surfaces. These templates are waiting for a big feature that we are adding to the query engine, also currently in testing.

New query functions:

TileInterpolateKriging - takes a drawing and the interpolation parameters and returns an interpolation object that can be used to produce tiles using kriging.

The heights are taken from Z values in the drawing geoms. Geoms without Z values are ignored. All geoms are converted to coordinates. Duplicate XY coordinates are ignored: if duplicate XY coordinates have different Z values, the function uses one of these values and ignores the other ones. (We can generate an error here, if necessary.)

The interpolation parameters are: radius for neighbor searches (kriging interpolates heights from closest coordinates), number of neighbors to use, interpolation model. If both the radius and the number of neighbors are zero or negative, kriging uses Voronoi neighbors (however many of them there would be). If both the radius and the number of neighbors are positive, kriging uses no more than the specified number of neighbors in the specified radius. (The behavior for the remaining two combinations is currently being adjusted.) The interpolation model can be set to a zero or negative value to autoselect the model, or to a specific value: 1 = linear, 2 = circular, 3 = exponential, 4 = gaussian, 5 = power, 6 = rational, 7 = spherical.

TileInterpolateKrigingPar - a parallel variant of TileInterpolateKriging, takes an additional parameter with the thread configuration.

TileInterpolateKrigingMedianPolish - takes the same parameters as TileInterpolateKriging and returns an interpolation object that can be used to produce tiles using kriging with median-polish (additional processing aiming to improve the quality of the interpolation).

TileInterpolateKrigingMedianPolishPar - a parallel variant of TileInterpolateKrigingMedianPolish, takes an additional parameter with the thread configuration.

TileInterpolateTriangulation - takes a drawing and returns an interpolation object that can be used to produce tiles using triangulation. (Like with kriging, all geoms are currently converted to coordinates. This is going to be changed, we will use constrained triangulation for lines and either convert areas to boundaries or ignore them.)

TileInterpolateTriangulationPar - a parallel variant of TileInterpolateTriangulation, takes an additional parameter with the thread configuration.

TileInterpolate - takes an interpolation object created by one of the above functions, a rectangle, and creates a tile with interpolated pixel values.

New service query functions:

CoordConvertRect - takes a coordinate converter object, a rect value, and the optional number of intermediate divisions, and converts the rect to the desired coordinate system. If the number of intermediate divisions is positive, the function creates a grid of intermediate coordinates covering the rect, projects all these coordinates and uses the results to compute the shape of the resulting rect more accurately (useful when the projection is curvilinear).

ValueSequenceTileXY - takes a rect value for the image, a tile size value, and a contained / touching boolean switch, and returns a table of XY indexes for tiles either completely within or with any part within the specified rect. The returned table includes the following fields: X - X coordinate of a tile, can be negative, Y - Y coordinate of a tile, can be negative, Rect - rect of a tile (it can be computed from X, Y and the tile size, but is provided regardless for convenience).

(The post after this one provides an example MXB with the query that illustrates how to use the above functions together.)


Matching brackets in text windows using Ctrl-] matches #...# delimiters used for dates in SQL.

The list of Standard coordinate systems no longer includes generic templates for coordinate systems like Orthographic which do not cover the entire Earth and need customization. (These templates come from Manifold 8 which used them as a starting point for creating a custom coordinate system. Leaving these templates in the list in the Standard tab was inviting the user to make a mistake of using them without further customization, so we removed them. All components set to use these templates will still continue to work / import / export as they should, the change affects only the list of displayed systems.)

The toolbar for a layout window or a map window collapses cursor modes into a single button with a dropdown menu. The icon on the button is set to the current cursor mode, the tooltip for the button displays the description of the current cursor mode.

The toolbar for a map window includes a button for locations, with a dropdown menu.

Switching the cursor mode in a map window cancels all pending editing. (For example, starting to enter coordinates for an area, then switching to tracker removes all entered coordinates and clears the Record pane.)

The toolbar no longer includes a button for the Help - About command.

The toolbar for a layout window no longer includes a button for the File - Print command (we have Ctrl-P to invoke this command, plus the change makes the toolbar line up with the toolbar for a map window). The File - Print Preview command has an icon.

Right-clicking a location in the Project pane allows applying the location to the active window using the new 'View in Active Window' command.

Saving a virtual layout, location, map or query as a persistent component puts the new component into the current path for the target data source in the Project pane. (This handles all cases where the new component goes to data source A but the Project pane has a completely different data source B selected, etc.)

Reading a TIFF interprets 3-channel floating-point images as RGB instead of BGR, 4-channel floating-point images as RGBA instead of BGRA, and sets the default value range to [0, 1].

Exporting data to a TIFF writes additional data for the projection into the file.

(Fix) Exporting data to a GDB no longer fails if the table contains a field with a reserved name (typically, OBJECTID).

(Fix) Exporting data to a GDB no longer produces wrong results if one or more of the fields fail to add. (The new code skips fields that fail to add and writes correct data into the remaining fields.)

Exporting data to an ECW or a JPEG2K writes image offset and scale info into the image file in addition to the projection info in the accompanying file.

Reading a KML file reads GroundOverlay data.

(Fix) Migrating images from MAP files created by Manifold 8 no longer sometimes shifts them 32 pixels up unnecessarily.

End of list.

8,093 post(s)
#19-Jun-18 22:48



7,972 post(s)
#19-Jun-18 15:46

Example of using kriging:

See the attached MXB.

The file contains an image (german_alps) and contours created on that image (german_alps_contours). There is also a query that takes the contours and interpolates them back into an image (puts the result into german_alps_kriging).

Run the query, wait for it to complete (again, the performance is not final, it is going to be faster, the query isn't even using threads), then check out the interpolated image.

Here is the text of the query:



-- put height into geoms

UPDATE [german_alps_contours Drawing] SET [geom] = GeomSetZ([geom][Value]);


-- drop old tiles

DELETE FROM [german_alps_kriging];


-- add new tiles

INSERT INTO [german_alps_kriging] ([x][y][tile])

SELECT [xy].[x][xy].[y],

  CASTV(TileInterpolate( -- interpolating a tile

    CALL TileInterpolateKriging([german_alps_contours Drawing], -1, -1, -1),

      -- interpolation object

    [xy].[rect]-- rect of the tile being interpolated


FROM CALL ValueSequenceTileXY(

  CAST (ComponentProperty([german_alps_kriging]'Rect')

    AS FLOAT64X4), -- rect in pixels

  CAST (ComponentProperty([german_alps_kriging Tiles]'FieldTileSize.Tile')

    AS FLOAT64X2), -- tile size

  false -- tiles with any part inside of the rect

AS [xy];


-- update intermediate levels

EXECUTE CALL TileUpdatePyramids([german_alps_kriging]);

The query first puts Z values computed for the contours into the geoms (UPDATE).

The query then clears the interpolated image by deleting all its tiles (DELETE) and repopulates it (INSERT).

The call to ValueSequenceTileXY in the FROM section of the INSERT produces the X / Y / Rect values for tiles to interpolate using the rect of the image and the tile size (both taken from properties of relevant components, could have been entered verbatim).

The tiles are produced using a call to TileInterpolate in the SELECT list. The interpolation object is created using TileInterpolateKriging, all parameters are set to 'auto' = Voronoi neighbors, automatic model. (The current build logs actual used values in the log, we will either remove the logging or tidy it up.) The produced tiles are converted to FLOAT32 using CASTV, because that's the type used by the image.

After producing the tiles, the query creates intermediate levels for the interpolated image (EXECUTE).

Templates coming soon.


8,093 post(s)
#19-Jun-18 22:49

Some quick notes.


I think Kriging has just become far more useful than it has ever been in Manifold 8, mostly because of speed. And this is only partly optimized and only using one thread--but wow! It is fast.

Speed is not just about getting the same amount done in less time, or even getting more done in the same time.

Even more important is that a disruptive speed increase like this allows experimentation (problem-solving, creativity, new questions, new products...) that previously would have been out of the question.

If I can try something 3 different ways in 10 minutes, then why not, even if two ways turn out to be (perhaps helpful, perhaps stupid) mistakes. On the other hand, if just one way would take nearly an hour, then I will probably not consider this avenue at all (bearing in mind not only the raw cost, but the price of a simple mistake), and I will instead do something less ambitious, perhaps less promising, more pedestrian.

So this is absolutely great and it is only going to get better with more optimizations and parallelism.


Also, the result of Kriging in this case is beautiful. Greate to be able to compare with the original.


I was surprised how fast the first UPDATE statement runs, using GeomSetZ. That's efficient.


Would it be possible to have Triangulation adjusted for contours (a.k.a. constrained triangulation + DEST) ready to try out, at around the same time as standard constrained triangulation?

Selfishly, that would make the most significant difference to my workflows.

By the same token, I will be able to give it a more thorough workout than (probably) anyone anywhere.


I really like the new ValueSequenceTileXY() construct. It's elegant, familiar in the context of SQL9, and much easier to use/read than a lower-level alternative (i.e. a set of ordinary ValueSequence() functions, probably with hard-coded values).

The way this works in with CAST(<string-property> AS <numeric-vector>) here is also really cool. You wired that sort of easy conversion up a long time ago; in a case like this it really pays off.


It's nice that you can use an alias with table like the one produced by CALL ValueSequenceTileXY(). I didn't know we could do that--obvious now given an example. On the other hand it's not necessary in this case, and the code might be clearer without it.


7,972 post(s)
#20-Jun-18 14:25

Thanks. :-)

We are working to add DEST together with the constrained triangulation, yes (shortly after it).

Good catch with the alias for the result of ValueSequenceTileXY() not being needed. I agree the query is clearer without it.

8,093 post(s)
#21-Jun-18 01:48

Very happy about that second para Adam, thank you for answering.

MSI_M92 post(s)
#28-Jun-18 16:04

Has there been any consideration to adding semi-variogram models to the Kriging routines in release 9?


7,972 post(s)
#28-Jun-18 16:25

We currently support 8 models, see the model parameter. You can either select a specific model or let the system determine the best model to use automatically.

We can absolutely add more models, if that's what you are asking.

MSI_M92 post(s)
#28-Jun-18 17:24

I need to read though the release notes more carefully it would seem!

My impression is that selection of models would satisfy our needs. I believe

with release 8 our hydrogeology folks would often use other software because

they couldn't control the semi-variogram in 8.

This is great - thanks for clarifying.


1,630 post(s)
#28-Jun-18 20:21

I might be misinterpreting what you are trying to say, but it sounds like your hydrology folks would like to see an interactive representation of the semivariogram for the dataset being interpolated:

I guess that this would allow the user to choose the model and parameters interactively as required. If you check out the log pane, M9 logs the model and parameters which is a good start.

Landsystems Ltd ... Know your land |


3,072 post(s)
#28-Jun-18 20:56

You can create a variogram with SQL. I demonstrate how to do this here.


456 post(s)
#02-Jul-18 23:14


join image "Because my dad promised me" interstellar from Manifold: Time by Stephen Baxter. power Math destruction


1,630 post(s)
#24-Jun-18 21:49

I think Kriging has just become far more useful than it has ever been in Manifold 8, mostly because of speed. And this is only partly optimized and only using one thread--but wow! It is fast.

Just to echo this. I hadn't realized just how much faster it already is (likely as I discounted Kriging from M8 because it was so much slower that some of the other interpolators). Anyway the same data and parameters in M8 have thus far been churning for 18hrs and 7mins in 8 and the process is 22% complete. In M9 on the same machine it took:

Query: [Interpolate Contours] (16.057 sec)


Landsystems Ltd ... Know your land |

Manifold User Community Use Agreement Copyright (C) 2007-2017 Manifold Software Limited. All rights reserved.