Add imsobel() and imfilter() (lepto) functions. (#1648)

* New tests fig

* Add 2 tests

* Export the 2 new funs.

* helper_mat2img(), cmap must have 4 columns.

* Add imsobel() and imfilter() (lepto) functions.

src/GMT.jl

@@ -167,8 +167,8 @@ export
 	VSdisp, mad, info, kmeans, pca, mosaic, quadbounds, quadkey, geocoder, getprovider,
 
 	bwhitmiss, binarize, bwperim, bwskell, isodata, padarray, rgb2gray, rgb2lab, rgb2YCbCr, rgb2ycbcr, grid2img,
-	img2grid, grays2cube, grays2rgb, imcomplement, imcomplement!, imerode, imdilate, imopen, imclose, imsegment,
-	imtophat, imbothat, imhdome, imhmin, imhmax, immorphgrad, imrankfilter, strel, 
+	img2grid, grays2cube, grays2rgb, imclose, imcomplement, imcomplement!, imdilate, imerode, imfilter, imopen, imsegment,
+	imsobel, imtophat, imbothat, imhdome, imhmin, imhmax, immorphgrad, imrankfilter, strel, 
 
 	imfill, imreconstruct, fillsinks, fillsinks!,
 

src/lepto_funs.jl

@@ -284,7 +284,7 @@ Here, a hole is defined as an area of dark pixels surrounded by lighter pixels.
 ### Examples
 ```julia
 # Example from Matlab imfill
-I = gmtread("C:\\programs\\MATLAB\\R2024b\\toolbox\\images\\imdata\\coins.png");
+I = gmtread("TESTSDIR * "assets/coins.png");
 Ibw = binarize(I);
 BW2 = imfill(Ibw);
 ```
@@ -947,6 +947,88 @@ function imrankfilter(I::GMTimage; width::Int=3, height::Int=0, rank=0.5)::GMTim
 	helper_morph(I, c, round(Int, rank * 1000), "rankf", nothing)
 end
 
+# ---------------------------------------------------------------------------------------------------
+"""
+    J = imsobel(I::GMTimage{<:UInt8, 2}; direction::Int=2)::GMTimage
+
+Sobel edge detecting filter.
+
+Returns a binary image BW containing 1s where the function finds edges in the grayscale or binary image I and 0s elsewhere.
+
+### Args
+- `I::GMTimage`: Input grayscale image.
+
+### Kwargs
+- `direction::Int=2`: Orientation flag: `0` means detect horizontal edges, `1` means vertical and `2` means all edges.
+
+### Returns
+A new `GMTimage` grayscale image `I` with the edge detection, edges are brighter.
+
+### Example
+
+Let us vectorize the rice grains in the image "rice.png". The result is not perfect because
+the grains in the image's edge are not closed and therefore not filled. And those get poligonized
+twice (in and outside) making the red lines look thicker (but they are not, they are just doubled).
+
+```julia
+I = gmtread(TESTSDIR * "assets/rice.png");
+J = imsobel(I);
+BW = binarize(J);
+BW = bwskell(BW);	# Skeletize to get a better approximation of the shapes.
+BW = imfill(BW);	# Fill to avoid "double" vectorization (outside and inside grain outline)
+D = polygonize(BW);
+grdimage(I, figsize=5)
+grdimage!(J, figsize=5, xshift=5.05)
+grdimage!(BW, figsize=5, xshift=-5.05, yshift=-5.05)
+grdimage!(I, figsize=5, xshift=5.05, plot=(data=D, lc=:red), show=true)
+```
+"""
+function imsobel(I::GMTimage{<:UInt8, 2}; direction::Int=2)::GMTimage
+	@assert 0 <= direction <= 2 "'direction' must be 0 (horizontal), 1 (vertical) or 2 (all edges)"
+	helper_morph(I, direction, 0, "sobel", nothing)
+end
+
+# ---------------------------------------------------------------------------------------------------
+"""
+	J = imfilter(I::GMTimage, kernel::Matrix{<:Real}; normalize::Int=1)::GMTimage
+
+Generic convolution filter.
+
+### Args
+- `I::GMTimage`: Input image. This can be a RGB or grayscale image.
+- `kernel::Matrix{<:Real}`: The filter kernel MxN matrix.
+
+### Kwargs
+- `normalize::Int=1`: Normalize the filter to unit sum. This is the default, unless ``sum(kernel) == 0``,
+   case in which we change it to 0.
+
+### Returns
+A new `GMTimage` of the same type as `I` with the filtered image.
+
+### Example
+
+Apply an 3x3 Mean Removal filter to the image "moon.png".
+
+```julia
+I = gmtread(TESTSDIR * "assets/moon.png");
+J = imfilter(I, [-1 -1 -1; -1 9 -1; -1 -1 -1]);
+grdimage(I, figsize=5)
+grdimage!(J, figsize=5, xshift=5, show=true)
+```
+"""
+function imfilter(I::GMTimage, kernel::Matrix{<:Real}; normalize::Int=1)::GMTimage
+	@assert (eltype(I) == UInt8) "'imfilter' is only available for UInt8 images"
+	(sum(kernel) == 0) && (normalize = 0)
+	bpp = (size(I,3) >= 3) ? 32 : 8
+	ppixI = img2pix(I, bpp)
+	kel = filtkernel(kernel)
+	if (size(I,3) == 1)  ppix = pixConvolve(ppixI.ptr, Ref(kel), 8, normalize)	# 8 => requesting to always return a 8 bpp image
+	else                 ppix = pixConvolveRGB(ppixI.ptr, Ref(kel))
+	end
+	(ppix == C_NULL) && (@warn("Convolution filter operation failed"); return GMTimage())
+	pix2img(Sppix(ppix))
+end
+
 # =====================================================================================================
 function helper_morph(I, hsize, vsize, tipo, sel)
 	bpp = (eltype(I) == Bool || I.range[6] == 1) ? 1 : 8
@@ -1005,6 +1087,8 @@ function helper_morph(I, hsize, vsize, tipo, sel)
 	elseif (tipo == "skell")
 		type, conn = bituncat2(hsize)
 		ppix = pixThinConnected(pI, type, conn, vsize)
+	elseif (tipo == "sobel")
+		ppix = pixSobelEdgeFilter(pI, hsize)
 	end
 	(ppix == C_NULL) && (@warn("Operation '$(tipo)' failed"); return GMTimage())
 	_I = pix2img(Sppix(ppix))
@@ -1067,6 +1151,15 @@ function strel(name::String, par1::Int, par2::Int=0)::Sel
 	strel(out, name=name)
 end
 
+# ---------------------------------------------------------------------------------------------------
+function filtkernel(mat::Matrix{<:Real})::L_Kernel
+	sy, sx = size(mat)
+	cx, cy = floor.(Int32, (size(mat))./2)
+	_mat = Float32.(mat)
+	data = [pointer(_mat[i,:]) for i in 1:size(_mat,1)]
+	L_Kernel(sy, sx, cy, cx, pointer(data))
+end
+
 function Base.show(io::IO, ::MIME"text/plain", sel::Sel)::Nothing
 	println(io, "Name: ", unsafe_string(sel.name))
 	mat = Matrix{Int32}(undef, sel.sy, sel.sx)

src/libleptonica.jl

@@ -99,6 +99,27 @@ struct Sel
     name::Cstring
 end
 
+"""
+    L_Kernel
+
+Kernel
+
+| Field | Note                            |
+| :---- | :------------------------------ |
+| sy    | kernel height                   |
+| sx    | kernel width                    |
+| cy    | y location of kernel origin     |
+| cx    | x location of kernel origin     |
+| data  | data[i][j] in [row][col] order  |
+"""
+struct L_Kernel
+	sy::Cint
+	sx::Cint
+	cy::Cint
+	cx::Cint
+	data::Ptr{Ptr{Cfloat}}
+end
+
 """
 When the garbage collector collects this object the associated Sppix object will be freed in the C library.
 """
@@ -202,3 +223,18 @@ function pixColorSegment(pixs, maxdist, maxcolors, selsize, finalcolors, debugfl
 end
 
 pixRankFilter(pixs, wf, hf, rank) = ccall((:pixRankFilter, liblept), Ptr{Pix}, (Ptr{Pix}, Cint, Cint, Cfloat), pixs, wf, hf, rank)
+pixSobelEdgeFilter(pixs, orientflag) = ccall((:pixSobelEdgeFilter, liblept), Ptr{Pix}, (Ptr{Pix}, Cint), pixs, orientflag)
+pixTwoSidedEdgeFilter(pixs, orientflag) = ccall((:pixTwoSidedEdgeFilter, liblept), Ptr{Pix}, (Ptr{Pix}, Cint), pixs, orientflag)
+
+function pixConvolve(pixs, kel, outdepth, normflag)
+	ccall((:pixConvolve, liblept), Ptr{Pix}, (Ptr{Pix}, Ptr{L_Kernel}, Cint, Cint), pixs, kel, outdepth, normflag)
+end
+function pixConvolveSep(pixs, kelx, kely, outdepth, normflag)
+    ccall((:pixConvolveSep, liblept), Ptr{Pix}, (Ptr{Pix}, Ptr{L_Kernel}, Ptr{L_Kernel}, Cint, Cint), pixs, kelx, kely, outdepth, normflag)
+end
+
+pixConvolveRGB(pixs, kel) = ccall((:pixConvolveRGB, liblept), Ptr{Pix}, (Ptr{Pix}, Ptr{L_Kernel}), pixs, kel)
+function pixConvolveRGBSep(pixs, kelx, kely)
+	ccall((:pixConvolveRGBSep, liblept), Ptr{Pix}, (Ptr{Pix}, Ptr{L_Kernel}, Ptr{L_Kernel}), pixs, kelx, kely)
+end
+

src/utils_types.jl

@@ -1193,9 +1193,10 @@ function helper_mat2img(mat, x::Vector{Float64}, y::Vector{Float64}, v::Vector{F
 	else
 		(size(mat,3) == 1) && (color_interp = "Gray")
 		if (!isempty(hdr) && (hdr[5] == 0 && hdr[6] == 1))	# A mask. Let's create a colormap for it
-			colormap = zeros(Int32, 256 * 3)
+			colormap = zeros(Int32, 256 * 4)
 			n_colors = 256;					# Because for GDAL we always send 256 even if they are not all filled
 			colormap[2] = colormap[258] = colormap[514] = 255
+			colormap[256*3+1:end] .= 255
 		else
 			colormap = zeros(Int32,3)		# Because we need an array
 		end

test/test_lepto_funs.jl

@@ -105,9 +105,13 @@ J = bwperim(I);
 I = gmtread(TESTSDIR * "assets/rice.png");
 J = imtophat(I, hsize=11, vsize=11);
 J = imbothat(I, hsize=11, vsize=11);
+J = imsobel(I);
 
 I = gmtread(TESTSDIR * "assets/bone.png");
 J = bwskell(I);
 
 I = gmtread(TESTSDIR * "assets/bunny_cenora.jpg");
 J = imsegment(I, colors=7);
+
+I = gmtread(TESTSDIR * "assets/moon.png");
+J = imfilter(I, [-1 -1 -1; -1 9 -1; -1 -1 -1]);