Subscribe to this thread
Home - General / All posts - A three band, RGB version of the 5x5 Unsharp filter
Dimitri


5,359 post(s)
#11-Feb-19 17:46

As much fun as it is to sharpen grayscale images, everybody really wants to sharpen RGB images. Here's a query that does that using the 5x5 Unsharp filter I posted in this thread. The RGB query uses two functions, the processtile function used in the earlier query, plus a new processRGB function that applies processtile to each of the channels (extracted using TileChannel) in turn and concatenating them back together into a three channel tile. This is a simple use of filters, using the same filter for each color channel and not weighting the channels differently when reassembling them.

Note that the CREATE TABLE statement has been adjusted to use FLOAT64x3 and float64x3 instead of just plain float64, since now we are using tiles that are three channel.

Tomorrow I'll show how to do a two filter Sobel or Prewitt filter, which the use of the two functions, compartmentalizing what's actually done to a tile within the processtile function, makes easy. I'll also post a project that can be downloaded.

The query text follows at the end of this post. But first, some visuals (click on the pictures to make them bigger):

Above is the original St Peters image in RGB.

Above is the image after unsharp sharpening.

Above is a big image of the Great Pyramid at Giza.

Above is the same image after unsharp sharpening.

The query text:

-- $manifold$

--

-- Template for custom RGB 3x3, 5x5, etc filters

-- using parameters for filter and radius, and two functions

--

-- RGB Unsharp 5x5

VALUE @filter TILE = StringJsonTile('[

-0.00391, -0.01563, -0.02344, -0.01563, -0.00391,

-0.01563, -0.06250, -0.09375, -0.06250, -0.01563,

-0.02344, -0.09375,  1.85938, -0.09375, -0.02344,

-0.01563, -0.06250, -0.09375, -0.06250, -0.01563,

-0.00391, -0.01563, -0.02344, -0.01563, -0.00391

]', 5, 5, 1, true);

VALUE @radius UINT8 = 2;

CREATE TABLE [RGB_Unsharp] (

  [X] INT32,

  [Y] INT32,

  [Tile] TILE,

  [mfd_id] INT64,

  INDEX [mfd_id_x] BTREE ([mfd_id]),

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

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

  PROPERTY 'FieldTileType.Tile' 'float64x3'

);

-- Do not forget to change the Rect values when using other images

CREATE IMAGE [RGB_Unsharp Image] (

  PROPERTY 'Table' '[RGB_Unsharp]',

  PROPERTY 'FieldTile' 'Tile',

  PROPERTY 'FieldX' 'X',

  PROPERTY 'FieldY' 'Y',

  PROPERTY 'Rect' '[ 0, 0, 1214, 862 ]'

);

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

  TileFilter(@t, @r, @f)

END;

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

TileChannelsConcat(

  processtile(TileChannel(@t,0), @r, @f),

TileChannelsConcat(

  processtile(TileChannel(@t,1), @r, @f),

  processtile(TileChannel(@t,2), @r, @f)))

END;

PRAGMA ('progress.percentnext' = '100');

INSERT INTO [RGB_Unsharp] (

  [X][Y],

  [Tile]

SELECT

  [X][Y],

  CASTV ((TileRemoveBorder(processRGB(TileCutBorder([st_peters], VectorMakeX2([X][Y]), @radius), @radius, @filter), @radius)) AS FLOAT64)

FROM [st_peters]

THREADS SystemCpuCount();

TABLE CALL TileUpdatePyramids([RGB_Unsharp Image]);

Enjoy!

Attachments:
rgb_unsharp_01.png
rgb_unsharp_02.png
rgb_unsharp_pyramids_01.png
rgb_unsharp_pyramids_02.png

drtees57 post(s)
#11-Feb-19 20:41

Wow! Impressive improvement. I noticed that the unsharp filter does add color artifacts to some of the distinct features on St. Peter's Basilica. Those artifacts would make the image less usable when zoomed in. However, for large scale images, the results are stunning.

Dimitri


5,359 post(s)
#12-Feb-19 07:56

That filter came from Taylor Petrick's site. I've used the example from that page in the illustrations below. (Here is the home page. He has a four part series on convolution filters.)

The example RGB unsharp filter is a simple example, applying the same 5x5 convolution matrix to all three channels the same way, with no post-processing afterward. That keeps the example simple enough so people can focus on the basics of using convolution filters with three channels.

Sites that allow you to apply a filter don't always tell you how they thereafter post-process the result, such as stretching or autocontrasting into a 0 to 255 output range. The example doesn't do any of that, to keep the query as simple as possible.

Artifacts occur when sharpening effects cause channels to get out of proportion with each other. That can be restored by applying full range in Style. You then need to apply medium autocontrast. (This particular image needs slightly less than medium, but that won't be an option for a point-and-click dialog until later this year when Manifold goes through a big image cycle and introduces numerous image editing / Photoshop-like controls.) If you do a Full Range and then a medium Autocontrast to the unsharpened St Peters image you eliminate the artifacts.

Here is an example using the example query with Petrick's sample image:

The unmodified image.

The final result of unsharp after applying Style controls to the channels to first do Full Range and then medium Autocontrast. It's slightly more contrasty than the 5x5 unsplash result in Petrick's page.

What you see immediately after applying the unsharp filter, before tinkering with Style.

The above is after applying Full Range. No more artifacts, but the image needs contrast adjusted to be brighter and more contrasty.

Attachments:
sample_unsharp_01.png
sample_unsharp_02.png
sample_unsharp_03.png
sample_unsharp_04.png

Dimitri


5,359 post(s)
#12-Feb-19 08:08

Here is the St Peters image after unsharp, followed by using Style to apply medium Autocontrast. In the other post, I applied Full Range first, to illustrate how the color artifacts arise from out-of-balance scaling, but you don't need to do that. You can simply apply Autocontrast and go straight to that.

The first image below uses medium Autocontrast, and the second uses 1.5 stdev autocontrast.

Attachments:
st_peters_postprocess.png
st_peters_postprocess_1.png

Dimitri


5,359 post(s)
#12-Feb-19 09:52

The project has been published in the product downloads page, as RGB_Filter_Examples.mxb

StanNWT
137 post(s)
#12-Feb-19 22:22

Hi Dimitri,

I used your 5x5 filter on an RGB satellite image which is in a Lambert Conformal Conic projection.

The "prj" file has the following information for the TIFF file:

PROJCS["NAD_1983_CSRS_Northwest_Territories_Lambert",

GEOGCS["GCS_North_American_1983_CSRS",DATUM["D_North_American_1983_CSRS",

SPHEROID["GRS_1980",6378137.0,298.257222101]],PRIMEM["Greenwich",0.0],

UNIT["Degree",0.0174532925199433]],PROJECTION["Lambert_Conformal_Conic"],

PARAMETER["False_Easting",0.0],PARAMETER["False_Northing",0.0],

PARAMETER["Central_Meridian",-112.0],PARAMETER["Standard_Parallel_1",62.0],

PARAMETER["Standard_Parallel_2",70.0],PARAMETER["Latitude_Of_Origin",0.0],

UNIT["Meter",1.0],AUTHORITY["EPSG",3581]]

The "TFW" file has the following information:

1.500000000000000

0.000000000000000

0.000000000000000

-1.500000000000000

-890657.2500000000

9225632.2500000000

When I run your routine on the data set, it comes out as psuedo mercator in the wrong place in the world, when I change the projection to lambert conformal conic it still is not in the right place in the world.

There was a series of GDAL errors when I imported that TIFF file and the rest in group of images.

2019-02-12 13:40:39 ** ASCII value for tag "GDALNoDataValue" contains null byte in value; value incorrectly truncated during reading due to implementation limitations

2019-02-12 13:40:39 ** ASCII value for tag "GDALNoDataValue" contains null byte in value; value incorrectly truncated during reading due to implementation limitations

2019-02-12 13:40:39 ** ASCII value for tag "GDALNoDataValue" contains null byte in value; value incorrectly truncated during reading due to implementation limitations

2019-02-12 13:40:39 ** ASCII value for tag "GDALNoDataValue" contains null byte in value; value incorrectly truncated during reading due to implementation limitations

2019-02-12 13:41:12 -- Import: S:\SPOT_150cm_NWT_Mosaic_2017\107B\107B07_rgb.tif (32.812 sec)

The nice thing is that the 2.57 GB RGB image only took ~92-94 seconds with different runs to process with this filter. I'm hoping to be able to use these and other filters at this speed to process imagery if need be and obviously I'd like the images to show up in the correct locations.

I'm including my edits of your code so you can see if there's any errors:

-- $manifold$

--

-- Template for custom RGB 3x3, 5x5, etc filters

-- using parameters for filter and radius, and two functions

--

-- RGB Unsharp 5x5

VALUE @filter TILE = StringJsonTile('[

-0.00391, -0.01563, -0.02344, -0.01563, -0.00391,

-0.01563, -0.06250, -0.09375, -0.06250, -0.01563,

-0.02344, -0.09375, 1.85938, -0.09375, -0.02344,

-0.01563, -0.06250, -0.09375, -0.06250, -0.01563,

-0.00391, -0.01563, -0.02344, -0.01563, -0.00391

]', 5, 5, 1, true);

VALUE @radius UINT8 = 2;

CREATE TABLE [107B07_RGB_Unsharp] (

[X] INT32,

[Y] INT32,

[Tile] TILE,

[mfd_id] INT64,

INDEX [mfd_id_x] BTREE ([mfd_id]),

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

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

PROPERTY 'FieldTileType.Tile' 'float64x3'

);

-- Do not forget to change the Rect values when using other images

CREATE IMAGE [107B07_RGB_Unsharp Image] (

PROPERTY 'Table' '[107B07_RGB_Unsharp]',

PROPERTY 'FieldTile' 'Tile',

PROPERTY 'FieldX' 'X',

PROPERTY 'FieldY' 'Y',

PROPERTY 'Rect' '[ 0, 0, 32868, 27993 ]'

);

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

TileFilter(@t, @r, @f)

) END;

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

TileChannelsConcat(

processtile(TileChannel(@t,0), @r, @f),

TileChannelsConcat(

processtile(TileChannel(@t,1), @r, @f),

processtile(TileChannel(@t,2), @r, @f)))

) END;

PRAGMA ('progress.percentnext' = '100');

INSERT INTO [107B07_RGB_Unsharp] (

[X], [Y],

[Tile]

) SELECT

[X], [Y],

CASTV ((TileRemoveBorder(processRGB(TileCutBorder([107B07_rgb], VectorMakeX2([X], [Y]), @radius), @radius, @filter), @radius)) AS FLOAT64)

FROM [107B07_rgb]

THREADS SystemCpuCount();

TABLE CALL TileUpdatePyramids([107B07_RGB_Unsharp Image]);

Dimitri


5,359 post(s)
#13-Feb-19 08:49

That's a lot to unpack.

I don't use GDAL all that much so I can't help you with that, although hearing about a series of GDAL errors when you used it to import the TIFF file does not give me confidence you imported the TIFF file into Manifold without errors. :-)

Split this into two sub-tasks:

1. Import the image correctly. If this doesn't get done right, all that follows will be done wrong.

a. Import the TIFF using Manifold.

b. Open it and verify the projection for the image is correct.

i. Check what the Component panel of the Content pane reports.

ii. create a map. Add a Bing Streets image server layer to the map.

iii. add the imported image to the map. Is it in the right place?

2. Adapt the filter to use your image.

a. The filter as written doesn't set the projection. You have to use Assign coordinate system to do that, if you want to do it manually.

b. If you want to do it automatically, Edit Query is your friend. :-)

i. open the source image.

ii. choose a transform, like Channels.

iii. Set the action button to Add Component

iv. Press Edit Query to see what it writes.

What you'll see is a query that includes in the CREATE TABLE statement the property that is required to create the new component with the correct coordinate system assigned. It will be something like this:

CREATE TABLE [TerrainDEM Curvature, Profile] (

  [X] INT32,

  [Y] INT32,

  [Tile] TILE,

  [mfd_id] INT64,

  INDEX [mfd_id_x] BTREE ([mfd_id]),

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

  PROPERTY 'FieldCoordSystem.Tile' '{ "CenterLat": 0, "CenterLon": 180,

 "Eccentricity": 0, "FalseEasting": 0, "FalseNorthing": 0, "LocalOffsetX":

 -3733200, "LocalOffsetY": -1866600, "LocalScaleX": 300, "LocalScaleY": 300,

 "MajorAxis": 1188300, "System""Cylindrical Equidistant""Unit""Meter",

 "UnitScale": 1, "UnitShort""m" }',

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

  PROPERTY 'FieldTileType.Tile' 'float64'

);

In the above, I had opened an image called TerrainDEM that was in cylindrical equidistant projection, I chose Curvature, Profile, and then I pressed Edit Query. The above is the CREATE TABLE part of that query. Note that it adds a FieldCoordSystem.Tile property. Since Manifold wrote it based on what the actual coordinate system of the tile is, we know it is correct.

Based on what Edit Query writes, you can adapt the query I posted by using the CREATE TABLE part of what Edit Query creates (adjusting the name of the table created, of course, to match what you use in the rest of the query).

Hope that helps!

Dimitri


5,359 post(s)
#13-Feb-19 11:09

It's poor practice to use absurdly large expressions (if the expression is big, better to use a query), but one way to avoid doing infrastructure is to have the Transform panel do the work for you: use an expression in the Expression tab of the Transform panel. I've posted an example in this new thread.

adamw


8,447 post(s)
#19-Feb-19 07:22

A slightly simpler way to set the projection than re-citing it in full - since the target image is going to use the same projection as the original image, we can just carry it over:

--SQL9

ALTER TABLE [107B07_RGB_Unsharp] (

  ADD PROPERTY 'FieldCoordSystem.Tile' ComponentCoordSystem([107B07_rgb])

);

Hope this helps.

Also, we are likely going to add unsharp as a built-in transform. Same for several other things.

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