From 8edf28e9b70aa81a58a414fff70f93c02ded12dc Mon Sep 17 00:00:00 2001 From: Joaquim Date: Sat, 1 Nov 2025 19:50:36 +0000 Subject: [PATCH] Fix several issues (but not all) with bwconncomp and image layouts. --- src/imgmorph/cc2bw.jl | 10 ++++-- src/lepto_funs.jl | 67 +++++++++++++++++++++++++++++------------ test/test_lepto_funs.jl | 19 +++++++----- 3 files changed, 65 insertions(+), 31 deletions(-) diff --git a/src/imgmorph/cc2bw.jl b/src/imgmorph/cc2bw.jl index 21b70daa8..25439eab1 100644 --- a/src/imgmorph/cc2bw.jl +++ b/src/imgmorph/cc2bw.jl @@ -1,13 +1,17 @@ """ + BW = cc2bw(cc::GMTConComp; obj2keep::Union{Int, Vector{Int}}=0) + +Convert connected components to binary image """ -function cc2bw(cc::GMTConComp) +function cc2bw(cc::GMTConComp; obj2keep::Union{Int, Vector{Int}}=0) # Create a binary image of the same size as the original I = mat2img(zeros(Bool, cc.image_size), x=cc.x, y=cc.y, inc=cc.inc, layout=cc.layout, is_transposed=(cc.layout[2] == 'R'), proj4=cc.proj4, wkt=cc.wkt, epsg=cc.epsg) # Mark pixels belonging to any connected component as true - for pixel_list in cc.pixel_list - for rc in pixel_list + plist = (obj2keep == 0) ? (1:cc.num_objects) : isvector(obj2keep) ? obj2keep : (obj2keep:obj2keep) + for k = 1:numel(plist) + for rc in cc.pixel_list[plist[k]] I[rc] = true end end diff --git a/src/lepto_funs.jl b/src/lepto_funs.jl index 547b9378c..367b112a7 100644 --- a/src/lepto_funs.jl +++ b/src/lepto_funs.jl @@ -1190,8 +1190,10 @@ Base.@kwdef struct GMTConComp proj4::String="" wkt::String="" epsg::Int=0 - bboxs::Vector{GMTdataset{Float64,2}}=GMTdataset{Float64,2}[] # The bounding boxes as datasets - pixel_list::Vector{Vector{Int32}}=Vector{Int32}[] # The list of pixel indices for each component + bbox::Vector{GMTdataset{Float64,2}}=GMTdataset{Float64,2}[] # The bounding boxes as datasets + pixel_list::Vector{Vector{Int32}}=Vector{Int32}[] # The list of pixel indices for each component + centroid::Matrix{Float64}=Matrix{Float64}(undef, 0, 0) # The centroids of each component + area::Vector{Float64}=Float64[] # The areas of each component end """ @@ -1220,8 +1222,10 @@ A `GMTConComp` structure with the following fields: - `proj4`: Projection definition of the image used to find the components. - `wkt`: Well-known text definition of the image used to find the components. - `epsg`: EPSG code of the image used to find the components. -- `bboxs`: The bounding boxes as a vector of GMTdataset. -- `pixel_list`: A vevtor of vectors with the list of linear pixel indices for each component. +- `bbox`: The bounding boxes as a vector of GMTdataset. +- `pixel_list`: A vector of vectors with the list of linear pixel indices for each component. +- `centroid`: A vector of Float64 Tuples with the x,y coordinates of the centroids for each component. +- `area`: A vector of Float64 with the areas of each component. This area is approximated by the mean lat for geog images. """ bwconncomp(mat::BitMatrix; conn::Int=8) = bwconncomp(mat2img(collect(mat)); conn=conn) bwconncomp(mat::Matrix{<:Bool}; conn::Int=8) = bwconncomp(mat2img(mat); conn=conn) @@ -1241,30 +1245,53 @@ function bwconncomp(I::GMTimage; conn::Int=8) #ppix = unsafe_load(unsafe_load(pixa).pix, 1) #return pix2img(Sppix(ppix)) - boxa = unsafe_load(pboxa) - D = Vector{GMTdataset{Float64,2}}(undef, boxa.n) + is_geog = isgeog(I) + cellsize = I.inc[1] * I.inc[2] + boxa = unsafe_load(pboxa) + D = Vector{GMTdataset{Float64,2}}(undef, boxa.n) + centroid = Matrix{Float64}(undef, boxa.n, 2) + area = Vector{Float64}(undef, boxa.n) pixel_list = [Int32[] for _ in 1:boxa.n] width, height = getsize(I) + is_col_maj = (I.layout[2] == 'C') + is_top_down = (I.layout[1] == 'T') + @inbounds for k = 1:boxa.n - b = unsafe_load(unsafe_load(boxa.box, k)) # A Box, x,y counting seem to be 0-based - y2 = I.y[Int32(width - b.y)] - y1 = I.y[Int32(width - b.y - b.h)] - x1, x2 = I.x[b.x + I.registration], I.x[b.x + b.w + I.registration] # Still not very clear where/why +1's are needed - D[k] = mat2ds([x1 y1; x1 y2; x2 y2; x2 y1; x1 y1]) # The bounding box - D[k].bbox = [x1, x2, y1, y2] # (xmin, xmax, ymin, ymax) - if (I.layout[2] == 'C') # Column major - jj = (I.layout[1] == 'T') ? (b.y+1:b.y+b.h) : (height-b.y:-1:height - b.y - b.h + 1) - @inbounds for i = b.x+1:b.x+b.w, j = jj - (I[j,i] != 0) && append!(pixel_list[k], (i-1) * height + j) + b = unsafe_load(unsafe_load(boxa.box, k)) # A Box, x,y counting seem to be 0-based + y2 = I.y[Int32(height - b.y) + 1] + y1 = I.y[Int32(height - b.y - b.h) + 1] + x1, x2 = I.x[b.x + I.registration], I.x[b.x + b.w + I.registration] + D[k] = mat2ds([x1 y1; x1 y2; x2 y2; x2 y1; x1 y1]) # The bounding box + D[k].bbox = [x1, x2, y1, y2] # (xmin, xmax, ymin, ymax) + half_pix_x, half_pix_y = I.registration * I.inc[1]/2, I.registration * I.inc[2]/2 + im, jm = 0.0, 0.0 + + if (is_col_maj) # Column major + jj = (is_top_down) ? (b.y+1:b.y+b.h) : (height-b.y:-1:height - b.y - b.h + 1) + @inbounds for i in b.x+1:b.x+b.w # Xaxis + @inbounds for j in jj + (I[j,i] != 0) && (append!(pixel_list[k], (i-1) * height + j); jm += j; im += i) + end end - else # Row major - @inbounds for j = b.y+1:b.y+b.h, i = b.x+1:b.x+b.w # Oddly here we do not need to care about 'T' or 'B' - (I[i,j] != 0) && append!(pixel_list[k], (j-1) * width + i) + n_px = length(pixel_list[k]) + t = (is_top_down) ? (I.range[4] - half_pix_y) - (jm / n_px - 1)*I.inc[2] : (I.range[3] + half_pix_x) + (jm / n_px - 1)*I.inc[2] + centroid[k,1:2] .= (I.range[1] + half_pix_x) + (im / n_px - 1)*I.inc[1], t + else # Row major + @inbounds for j = b.y+1:b.y+b.h # Yaxis. + @inbounds for i = b.x+1:b.x+b.w + (I[i,j] != 0) && (append!(pixel_list[k], (j-1) * width + i); im += i; jm += j) + end end + n_px = length(pixel_list[k]) + t = (is_top_down) ? (I.range[2] - half_pix_y) - (jm / n_px - 1)*I.inc[1] : (I.range[1] + half_pix_y) + (jm / n_px - 1)*I.inc[1] + centroid[k,1:2] .= (I.range[3] + half_pix_x) + (im / n_px - 1)*I.inc[2], t end + lat_fact = (is_geog) ? cosd(centroid[k][2]) : 1.0 # Big approximation for the geog case but exact solution is too complex here + area[k] = n_px * cellsize * lat_fact end set_dsBB!(D, true) # Add the bounding box - GMTConComp(conn, pcount[], size(I), I.range, I.inc, I.registration, I.x, I.y, I.layout, I.proj4, I.wkt, I.epsg, D, pixel_list) + GMTConComp(conn, pcount[], size(I), I.range, I.inc, I.registration, I.x, I.y, I.layout, I.proj4, I.wkt, I.epsg, + D, pixel_list, centroid, area) end function get_ppixI(I) diff --git a/test/test_lepto_funs.jl b/test/test_lepto_funs.jl index 885ad8761..0636e5daa 100644 --- a/test/test_lepto_funs.jl +++ b/test/test_lepto_funs.jl @@ -62,7 +62,7 @@ end @test filled_gmt.image[1, 1] == 255 end -#= +## @testset "bwconncomp tests" begin I = gdalread(TESTSDIR * "assets/packman.png"); cc1 = bwconncomp(I); @@ -70,13 +70,16 @@ end I = gdalread(TESTSDIR * "assets/packman.png", layout="BCBa"); cc2 = bwconncomp(I); @test cc2.num_objects == 2 - @test cc1.bboxs == cc2.bboxs - #I = gdalread(TESTSDIR * "assets/packman.png", layout="TCBa"); # Need to fix bug in gdalread for this to work - #cc2 = bwconncomp(I); - #@test cc2.num_objects == 2 - #@test cc1.bboxs == cc2.bboxs + @test cc1.bbox == cc2.bbox + I = gdalread(TESTSDIR * "assets/packman.png", layout="TCBa"); + cc2 = bwconncomp(I); + @test cc2.num_objects == 2 + @test cc1.bbox == cc2.bbox # This guy fails - mat = zeros(Bool,8,9); mat[2:5,2:3] .= true; mat[4:5,4] .= true; mat[5:8,6:8] .= true; mat[4,7] = true; + # LEPTONICA is BUGGED. This trims last column of second object + #mat = zeros(Bool,8,10); mat[2:5,2:3] .= true; mat[4:5,4] .= true; mat[5:8,7:9] .= true; mat[4,8] = true; + # But this one (moving second obj one pixel to the left) works + mat = zeros(Bool,8,10); mat[2:5,2:3] .= true; mat[4:5,4] .= true; mat[5:8,6:8] .= true; mat[4,7] = true; gdalwrite("c2.png", mat2img(mat)) I = gdalread("c2.png"); # Comes in with layout TRBa cc = bwconncomp(I); @@ -94,7 +97,7 @@ end I2 = cc2bw(cc); @test I == I2 end -=# +## G = peaks(); G2 = fillsinks(G);