SQL Example: Process Images using Dual 3x3 Filters

This topic is a continuation of the SQL Example: Process Images with 3x3 Filters  topic.  Please read that topic before continuing, as this topic picks up directly where that last topic finished.

 

The discussion in this topic is the second in a series of three topics, beginning with the SQL Example: Process Images with 3x3 Filters  topic and continuing after this topic in the  SQL Example: Process RGB Images using Matrix Filters topic.

 

 This topic uses the same images as examples, including  the st_peters_grayscale image shown below, slightly zoomed in.

 

More Sophisticated Processing

One reason we use a processtile function is to compartmentalize the processing of each tile within a function.   That makes it easier to do more extensive processing than to simply apply TileFilter once.

 

Consider our sample query from the end of the SQL Example: Process Images with 3x3 Filters  topic, as adjusted to use Sobel vertical edge filter with a 3x3 matrix and the original st_peters_grayscale image:

 

-- prepare begin

 

CREATE TABLE [Custom] (

  [X] INT32,

  [Y] INT32,

  [mfd_id] INT64,

  [Tile] TILE,

  INDEX [mfd_id_x] BTREE ([mfd_id]),

  INDEX [X_Y_Tile_x] RTREE ([X], [Y], [Tile] TILESIZE (128, 128) TILETYPE FLOAT64),

  PROPERTY 'FieldCoordSystem.Tile' ComponentFieldCoordSystem([st_peters_grayscale Tiles], 'Tile'),

  PROPERTY 'FieldTileSize.Tile' '[ 128, 128 ]',

  PROPERTY 'FieldTileType.Tile' 'float64'

);

CREATE IMAGE [Custom Image] (

  PROPERTY 'Table' '[Custom]',

  PROPERTY 'FieldX' 'X',

  PROPERTY 'FieldY' 'Y',

  PROPERTY 'FieldTile' 'Tile',

  PROPERTY 'Rect' CAST(ComponentFieldBounds([st_peters_grayscale Tiles], 'Tile', FALSE) AS NVARCHAR)

);

 

-- prepare end

 

VALUE @image TABLE = CALL ComponentFieldImage([st_peters_grayscale Tiles], 'Tile');

 

-- Vertical Edge Detect Single Sobel

VALUE @filterDef TILE = StringJsonTile('[ 1, 0, -1, 2, 0, -2, 1, 0, -1 ]', 3, 3, 1, true);

 

VALUE @radius UINT8 = 1;

 

FUNCTION processtile(@t TILE, @r UINT8, @f TILE) TILE AS (

  TileFilter(@t, @r, @f)

) END;

 

DELETE FROM [Custom];

INSERT INTO [Custom] (

  [X], [Y],

  [Tile]

) SELECT

  [X], [Y],

  CASTV ((TileRemoveBorder(processtile(

    TileCutBorder(@image, VectorMakeX2([X], [Y]), @radius), @radius, @filterDef

    ), @radius)) AS FLOAT64)

FROM [st_peters_grayscale Tiles] THREADS SystemCpuCount();

 

TABLE CALL TileUpdateFieldPyramids([Custom], 'Tile');

ALTER IMAGE [Custom Image] (

  ADD PROPERTY 'Rect' Coalesce(

    CAST(ComponentFieldBounds([Custom], 'Tile', TRUE) AS NVARCHAR),

    ComponentProperty([Custom Image], 'Rect'))

);

 

In the above query, the processtile function is simply:

 

FUNCTION processtile(@t TILE, @r UINT8, @f TILE) TILE AS (

  TileFilter(@t, @r, @f)

) END;

 

It calls TileFilter once to compute the resulting tile.   A more sophisticated approach to Sobel edge detection using convolution matrices is to use two different filter matrices, computing using one matrix to get a result, then computing using the other matrix to get a second result, and to then combine the two using the formula:

 

((first result)^2 + (second result)^2)^0.5

 

That is, we square the first result, add that to the square of the second result, and then we take the square root of the sum.   If we have two different filters to use, we can re-write the processtile function as:

 

FUNCTION processtile(@t TILE, @r UINT8, @f1 TILE, @f2 Tile) TILE AS (

  (TileFilter(@t, @r, @f1)^2 + TileFilter(@t, @r, @f2)^2)^0.5

) END;

 

Adjusting the sample query

In addition to changing the processtile function to take two filters as arguments and to do the above computation,, we will also have to re-write the CASTV SQL lines which call processtile, to pass not one but two filters to the function, and we will have to provide two VALUE statements for the two filters to use.   The sample query now looks like:

 

-- prepare begin

 

CREATE TABLE [Custom] (

  [X] INT32,

  [Y] INT32,

  [mfd_id] INT64,

  [Tile] TILE,

  INDEX [mfd_id_x] BTREE ([mfd_id]),

  INDEX [X_Y_Tile_x] RTREE ([X], [Y], [Tile] TILESIZE (128, 128) TILETYPE FLOAT64),

  PROPERTY 'FieldCoordSystem.Tile' ComponentFieldCoordSystem([st_peters_grayscale Tiles], 'Tile'),

  PROPERTY 'FieldTileSize.Tile' '[ 128, 128 ]',

  PROPERTY 'FieldTileType.Tile' 'float64'

);

CREATE IMAGE [Custom Image] (

  PROPERTY 'Table' '[Custom]',

  PROPERTY 'FieldX' 'X',

  PROPERTY 'FieldY' 'Y',

  PROPERTY 'FieldTile' 'Tile',

  PROPERTY 'Rect' CAST(ComponentFieldBounds([st_peters_grayscale Tiles], 'Tile', FALSE) AS NVARCHAR)

);

 

-- prepare end

 

VALUE @image TABLE = CALL ComponentFieldImage([st_peters_grayscale Tiles], 'Tile');

 

-- Sobel Edge Detect Filter 1

VALUE @filter1 TILE = StringJsonTile('[ -1, 0, 1, -2, 0, 2, -1, 0, 1 ]', 3, 3, 1, true);

 

-- Sobel Edge Detect Filter 2

VALUE @filter2 TILE = StringJsonTile('[ -1, -2, -1, 0, 0, 0, 1, 2, 1 ]', 3, 3, 1, true);

 

VALUE @radius UINT8 = 1;

 

FUNCTION processtile(@t TILE, @r UINT8, @f1 TILE, @f2 Tile) TILE AS (

  (TileFilter(@t, @r, @f1)^2 + TileFilter(@t, @r, @f2)^2)^0.5

) END;

 

DELETE FROM [Custom];

INSERT INTO [Custom] (

  [X], [Y],

  [Tile]

) SELECT

  [X], [Y],

  CASTV ((TileRemoveBorder(processtile(

    TileCutBorder(@source_image, VectorMakeX2([X], [Y]), @radius), @radius, @filter1, @filter2

    ), @radius)) AS FLOAT64)

FROM [st_peters_grayscale Tiles] THREADS SystemCpuCount();

 

TABLE CALL TileUpdateFieldPyramids([Custom], 'Tile');

ALTER IMAGE [Custom Image] (

  ADD PROPERTY 'Rect' Coalesce(

    CAST(ComponentFieldBounds([Custom], 'Tile', TRUE) AS NVARCHAR),

    ComponentProperty([Custom Image], 'Rect'))

);

 

The changes made have been indicated in blue color.   We now use two filter arguments, @filter1 and @filter2, in the expression that calls processtile, a simple change given the compartmentalization of function within the function.

 

 

The result shows the excellent detection of edges a two-filter Sobel edge detection process can accomplish.  

Dual Prewitt Filters

Prewitt filters also use two filters with the same "add the squares of the results and take the square root" formula used above, but they use different filter matrices.    We can replace the Sobel filters in our sample query with Prewitt filters to get:

 

 

-- prepare begin

 

CREATE TABLE [Custom] (

  [X] INT32,

  [Y] INT32,

  [mfd_id] INT64,

  [Tile] TILE,

  INDEX [mfd_id_x] BTREE ([mfd_id]),

  INDEX [X_Y_Tile_x] RTREE ([X], [Y], [Tile] TILESIZE (128, 128) TILETYPE FLOAT64),

  PROPERTY 'FieldCoordSystem.Tile' ComponentFieldCoordSystem([st_peters_grayscale Tiles], 'Tile'),

  PROPERTY 'FieldTileSize.Tile' '[ 128, 128 ]',

  PROPERTY 'FieldTileType.Tile' 'float64'

);

CREATE IMAGE [Custom Image] (

  PROPERTY 'Table' '[Custom]',

  PROPERTY 'FieldX' 'X',

  PROPERTY 'FieldY' 'Y',

  PROPERTY 'FieldTile' 'Tile',

  PROPERTY 'Rect' CAST(ComponentFieldBounds([st_peters_grayscale Tiles], 'Tile', FALSE) AS NVARCHAR)

);

 

-- prepare end

 

VALUE @image TABLE = CALL ComponentFieldImage([st_peters_grayscale Tiles], 'Tile');

 

-- Prewitt Edge Detect Filter 1

VALUE @filter1 TILE = StringJsonTile('[ -1, 0, 1, -1, 0, 1, -1, 0, 1 ]', 3, 3, 1, true);

 

-- Prewitt Edge Detect Filter 2

VALUE @filter2 TILE = StringJsonTile('[ -1, -1, -1, 0, 0, 0, 1, 1, 1 ]', 3, 3, 1, true);

 

VALUE @radius UINT8 = 1;

 

FUNCTION processtile(@t TILE, @r UINT8, @f1 TILE, @f2 Tile) TILE AS (

  (TileFilter(@t, @r, @f1)^2 + TileFilter(@t, @r, @f2)^2)^0.5

) END;

 

DELETE FROM [Custom];

INSERT INTO [Custom] (

  [X], [Y],

  [Tile]

) SELECT

  [X], [Y],

  CASTV ((TileRemoveBorder(processtile(

    TileCutBorder(@image, VectorMakeX2([X], [Y]), @radius), @radius, @filter1, @filter2

    ), @radius)) AS FLOAT64)

FROM [st_peters_grayscale Tiles] THREADS SystemCpuCount();

 

TABLE CALL TileUpdateFieldPyramids([Custom], 'Tile');

ALTER IMAGE [Custom Image] (

  ADD PROPERTY 'Rect' Coalesce(

    CAST(ComponentFieldBounds([Custom], 'Tile', TRUE) AS NVARCHAR),

    ComponentProperty([Custom Image], 'Rect'))

);

 

We run that query to see what difference using Prewitt filters makes:

 

 

Prewitt filters also do a fine job of edge detection, providing a slightly more contrasty result than Sobel filters.

 

Yet More Improvements

The next frontier for even more improvements is to take the sample queries we have created for single-channel images and to adapt them for operations on three-channel, RGB images.  The discussion continues in the SQL Example: Process RGB Images using Matrix Filters topic.

 

Videos

Manifold Viewer - How Matrix Filters Work - The easy, simple way to learn how filters work! Watch this action-packed video using Manifold Viewer that illustrates how matrix filters, also known as convolution filters, work. In addition to explaining filters, the video provides a real-life look at simple Manifold techniques for moving objects around in drawings using the Shift transform, and fast and easy use of Selection and tables to quickly put desired values into attributes. Sound technical? Sure, but in a very easy and simple way.

 

Manifold Viewer - Create Custom GPU Accelerated Filters in Seconds - A technical video using the free Viewer showing how to create your own, fully custom, fully GPU-parallel, convolution matrix filters, including Emboss, Sobel, Prewitt, and Kirsch edge detection and many more, for use in Viewer or Release 9. Modify the spatial SQL examples in the downloadable example project to specify a custom matrix and in seconds your custom filter can do image processing at GPU-parallel speeds. Viewer is read-only, but you can copy and paste the query text for custom filters to and from Notepad or any other text editor. Download the Custom_Filter_Examples.mxb sample project from the Examples page on the Manifold website to try out the video in Viewer or Release 9.

 

Manifold Viewer - Speed Demo with 1280 GPU Cores - 2 Minutes vs 5 Days - Watch the free Manifold Viewer run CPU parallel and GPU parallel with 8 CPU cores and 1280 GPU cores to absolutely crush a job, doing in 2 minutes what would take non-GPU-parallel software 5 days. The video shows Viewer instantly pop open a 4 GB project that contains a huge, multi-gigabyte terrain elevation surface for Crater Lake, Oregon. With a point and click - no parallel code required - we compute the mean curvature at each pixel of the surface using a 7x7 matrix in under two minutes. We combine that with the original surface for enhanced hill shaded effects to better see details. Using an 11x11 matrix takes just over two minutes, a huge computation that takes days in non-GPU-parallel GIS packages.

 

See Also

Images

 

Tables

 

Data Types

 

How Matrix Filters Work

 

Command Window

 

SQL Functions

 

SQL Example: Process Images with 3x3 Filters -  Shows a step-by-step example of developing an SQL query that takes a query written by the Edit Query button and then modifies that query into a general purpose query that can apply any 3x3 filter.   This makes it easy to use matrix filters we find on the web for custom image processing.   We extend the query by using parameters and adding a function, and then show how it can be adapted to use a 5x5 filter.

 

SQL Example: Process RGB Images using Matrix Filters - A continuation of the above topic and this topic, extending the example query to process three channel, RGB images.

 

SQL Example: Create NDVI Displays - How to create a query that creates an NDVI display from a four-band NAIP image, with tips and tricks on how to copy and paste existing information to get the result we want.

 

SQL Example: Create Topographic Position Index TPI Displays - In this example, we use a few short lines of SQL to create a Topographic Position Index (TPI) display.  TPI characterizes the undulations of a terrain elevation surface.  TPI value above zero show locations that are higher than then average of immediately surrounding terrain, and thus tend to show ridges.   TPI values below zero show locations that are lower than the average of immediately surrounding terrain, and thus tend to show valleys.  TPI values that are zero show areas of constant slope.

 

Example: Enhance Terrain with Curvatures -  We enhance a terrain showing Crater Lake, Oregon, by using mean curvature calculation to bring out details.   The example uses a 4 GB project containing a large terrain elevation surface.  Using a point-and-click dialog with no SQL, we apply automatic CPU parallelism and GPU parallelism to absolutely crush a task in two and a half minutes that would take non-parallel software days.

 

Example: Rearrange Channels using an Expression - We use a simple expression in the Transform pane to rearrange the order of channels within the data.