|
| 1 | +function find_histo_limits(In, thresholds=nothing, width=20, hst_::Matrix{Float64}=Matrix{Float64}(undef,0,2)) |
| 2 | + # Find the histogram limits of a UInt16 GMTimage that allow to better stretch the histogram |
| 3 | + # THRESHOLDS is an optional Tuple input containing the left and right histo thresholds, in percentage, |
| 4 | + # between which the histogram values will be retained. Defaults are (0.1, 0.4) percent. Note, this |
| 5 | + # assumes the histogram follows some sort of Gaussian distribution. It it's flat, shit occurs. |
| 6 | + # WIDTH is bin width used to obtain a rough histogram that is used to compute the limits. |
| 7 | + if (isa(In, Array{UInt16,3}) || isa(In, Array{UInt8,3})) |
| 8 | + L1 = find_histo_limits(view(In, :, :, 1), thresholds, width) |
| 9 | + L2 = find_histo_limits(view(In, :, :, 2), thresholds, width) |
| 10 | + L3 = find_histo_limits(view(In, :, :, 3), thresholds, width) |
| 11 | + return (L1[1], L1[2], L2[1], L2[2], L3[1], L3[2]) |
| 12 | + end |
| 13 | + hst = (isempty(hst_)) ? loc_histo(In, "", string(width), "")[1] : hst_ |
| 14 | + if (size(hst, 1) > 10) |
| 15 | + all(hst[2:5,2] .== 0) && (hst[1,2] = 0) # Here we always check for high counts in zero bin |
| 16 | + # Some processed bands leave garbage on the low DNs and that fouls our detecting algo. So check more |
| 17 | + ((hst[1,2] != 0) && hst[1,2] > 100 * mean(hst[2:10,2])) && (hst[1,2] = 0) # Ad-hoc numbers |
| 18 | + # Next is for the case of VIIRS bands that have nodata = 65535 but when called from 'truecolor' |
| 19 | + # it only passes a band view (the array) and we have access to image's header here. |
| 20 | + (width == 20 && eltype(In) == UInt16 && size(hst,1) == 3277) && (hst[end, 2] = 0) |
| 21 | + end |
| 22 | + max_ = maximum(hst, dims=1)[2] |
| 23 | + (max_ == 0) && error("This histogram had nothing but countings ONLY in first bin. No point to proceed.") |
| 24 | + thresh_l::Float64 = 0.001; thresh_r::Float64 = 0.004 |
| 25 | + if (isa(thresholds, Tuple) && length(thresholds) == 2) |
| 26 | + thresh_l, thresh_r = thresholds[:] ./ 100 |
| 27 | + end |
| 28 | + thresh_l *= max_ |
| 29 | + thresh_r *= max_ |
| 30 | + kl = 1; kr = size(hst, 1) |
| 31 | + while (hst[kl,2] == 0 || hst[kl,2] < thresh_l) kl += 1 end |
| 32 | + while (hst[kr,2] == 0 || hst[kr,2] < thresh_r) kr -= 1 end |
| 33 | + #return Int(hst[kl,1]), Int(min(hst[kr,1] + width, hst[end,1])) |
| 34 | + return hst[kl,1], min(hst[kr,1] + width, hst[end,1]) |
| 35 | +end |
| 36 | + |
| 37 | +# --------------------------------------------------------------------------------------------------- |
| 38 | +function loc_histo(in, cmd::String="", opt_T::String="", opt_Z::String="") |
| 39 | + # Very simple function to compute histograms of images (integers) |
| 40 | + # We put the countings in a Mx2 arrray to trick GMT (pshistogram) to think it's recieving a weighted input. |
| 41 | + (!isa(in[1], UInt16) && !isa(in[1], UInt8)) && error("Only UInt8 or UInt16 image types allowed here") |
| 42 | + |
| 43 | + inc::Float64 = (opt_T != "") ? parse(Float64, opt_T) : 1.0 |
| 44 | + (inc <= 0) && error("Bin width must be a number > 0 and no min/max") |
| 45 | + |
| 46 | + n_bins::Int = (isa(in[1], UInt8)) ? 256 : Int(ceil((maximum(in) + 1) / inc)) # For UInt8 use the full [0 255] range |
| 47 | + hst = zeros(n_bins, 2) |
| 48 | + pshst_wall!(in, hst, inc, n_bins) |
| 49 | + |
| 50 | + cmd = (opt_Z == "") ? cmd * " -Z0" : cmd * opt_Z |
| 51 | + (!occursin("+w", cmd)) && (cmd *= "+w") # Pretending to be weighted is crutial for the trick |
| 52 | + |
| 53 | + return hst, cmd * " -T0/$(n_bins * inc)/$inc" |
| 54 | +end |
| 55 | + |
| 56 | +# --------------------------------------------------------------------------------------------------- |
| 57 | +function pshst_wall!(in, hst, inc, n_bins::Int) |
| 58 | + # Function barrier for type instability. With the body of this in the calling fun the 'inc' var |
| 59 | + # introduces a mysterious type instability and execution times multiply by 3. |
| 60 | + if (inc == 1) |
| 61 | + @inbounds Threads.@threads for k = 1:numel(in) hst[in[k] + 1, 2] += 1 end |
| 62 | + else |
| 63 | + @inbounds Threads.@threads for k = 1:numel(in) hst[Int(floor(in[k] / inc) + 1), 2] += 1 end |
| 64 | + end |
| 65 | + (isa(in, GItype) && in.nodata == typemax(eltype(in))) && (hst[end] = 0) |
| 66 | + @inbounds Threads.@threads for k = 1:n_bins hst[k,1] = inc * (k - 1) end |
| 67 | + return nothing |
| 68 | +end |
| 69 | + |
| 70 | +# --------------------------------------------------------------------------------------------------- |
| 71 | +# This version computes the histogram for a UInt8 image band with a bin width of 1 |
| 72 | +histogray(img::GMTimage{<:UInt8}; band=1) = histogray(view(img.image, :, :, band)) |
| 73 | +function histogray(img::AbstractMatrix{UInt8}) |
| 74 | + edges, counts = 0:255, fill(0, 256) |
| 75 | + Threads.@threads for v in img |
| 76 | + @inbounds counts[v+1] += 1 |
| 77 | + end |
| 78 | + return counts, edges |
| 79 | +end |
| 80 | + |
| 81 | +# --------------------------------------------------------------------------------------------------- |
| 82 | +function hst_floats(arg1, opt_T::String=""; min_max=(0.0, 0.0)) |
| 83 | + # Compute the histogram of a grid or matrix |
| 84 | + # Made a separate function to let it be called from rescale() and thus avoid calling the main histogram() |
| 85 | + # that seems to be be too havy (at least according to JET) |
| 86 | + _min_max::Tuple{Float64,Float64} = (isa(arg1, GMTgrid)) ? |
| 87 | + (arg1.range[5], arg1.range[6]) : (min_max != (0.0, 0.0) ? min_max : Float64.(extrema_nan(arg1))) |
| 88 | + if (opt_T != "") |
| 89 | + inc = parse(Float64, opt_T) + eps() # + EPS to avoid the extra last bin at right with 1 count only |
| 90 | + n_bins = Int(ceil((_min_max[2] - _min_max[1]) / inc)) |
| 91 | + else |
| 92 | + n_bins = Int(ceil(sqrt(length(arg1)))) |
| 93 | + reg = isa(arg1, GMTgrid) ? (1 - arg1.registration) : 1 # When called from RemoteS arg1 is a view of a layer. |
| 94 | + inc = (_min_max[2] - _min_max[1]) / (n_bins - reg) + eps() |
| 95 | + end |
| 96 | + (!isa(inc, Real) || inc <= 0) && error("Bin width must be a > 0 number and no min/max") |
| 97 | + hst = zeros(n_bins, 2) |
| 98 | + have_nans = false |
| 99 | + if (eltype(arg1) <: AbstractFloat) # Float arrays can have NaNs |
| 100 | + have_nans = !(isa(arg1, GMTgrid) && arg1.hasnans == 1) |
| 101 | + have_nans && (have_nans = any(!isfinite, arg1)) |
| 102 | + end |
| 103 | + |
| 104 | + _inc = inc + 10eps() # To avoid cases when index computing fall of by 1 |
| 105 | + if (have_nans) # If we have NaNs in the grid, we need to take a slower branch |
| 106 | + @inbounds for k = 1:numel(arg1) |
| 107 | + !isnan(arg1[k]) && (hst[Int(floor((arg1[k] - _min_max[1]) / _inc) + 1), 2] += 1) |
| 108 | + end |
| 109 | + else |
| 110 | + @inbounds for k = 1:numel(arg1) hst[Int(floor((arg1[k] - _min_max[1]) / _inc) + 1), 2] += 1 end |
| 111 | + end |
| 112 | + @inbounds for k = 1:n_bins hst[k,1] = _min_max[1] + inc * (k - 1) end |
| 113 | + return hst, inc, _min_max |
| 114 | +end |
| 115 | + |
| 116 | +# --------------------------------------------------------------------------------------------------- |
| 117 | +function binmethod(d::Dict, cmd::String, X, is_datetime::Bool) |
| 118 | + # Compute bin width for a series of binning alghoritms or intervals when X (DateTime) comes in seconds |
| 119 | + val::String = ((val_ = find_in_dict(d, [:binmethod :BinMethod])[1]) !== nothing) ? lowercase(string(val_)) : "" |
| 120 | + min_max = (zero(eltype(X)), zero(eltype(X))) |
| 121 | + (!is_datetime) && (min_max = extrema(X)) # X should already be sorted but don't trust that |
| 122 | + if (val == "") |
| 123 | + if (!is_datetime) |
| 124 | + val = "sqrt" |
| 125 | + else |
| 126 | + min_max = extrema(X) # X should already be sorted but don't trust that |
| 127 | + rng = (min_max[2] - min_max[1]) |
| 128 | + if (rng < 150) val = "second" |
| 129 | + elseif (rng / (60) < 150) val = "minute" |
| 130 | + elseif (rng / (3600) < 150) val = "hour" |
| 131 | + elseif (rng / (86400) < 150) val = "day" |
| 132 | + elseif (rng / (86400 * 7) < 150) val = "week" |
| 133 | + elseif (rng / (86400 * 31) < 150) val = "month" |
| 134 | + else val = "year" |
| 135 | + end |
| 136 | + end |
| 137 | + end |
| 138 | + |
| 139 | + n_bins = 0.0; bin = 0 |
| 140 | + if (val == "scott") n_bins = 3.5 .* std(X) .* length(X)^(-1/3) |
| 141 | + elseif (val == "fd") n_bins = 2 .* IQR(X) .* length(X)^(-1/3) |
| 142 | + elseif (val == "sturges") n_bins = ceil.(1 .+ log2.(length(X))) |
| 143 | + elseif (val == "sqrt") n_bins = ceil.(sqrt(length(X))) |
| 144 | + elseif (val == "year") bin = 86400 * 365.25 |
| 145 | + elseif (val == "month") bin = 86400 * 31 |
| 146 | + elseif (val == "week") bin = 86400 * 7 |
| 147 | + elseif (val == "day") bin = 86400 |
| 148 | + elseif (val == "hour") bin = 3600 |
| 149 | + elseif (val == "minute") bin = 60 |
| 150 | + elseif (val == "second") bin = 1 |
| 151 | + elseif (!is_datetime) error("Unknown BinMethod $val") |
| 152 | + end |
| 153 | + if (bin == 0) |
| 154 | + bin = (min_max[2] - min_max[1]) / n_bins # Should be made a "pretty" number? |
| 155 | + end |
| 156 | + return @sprintf("%.12g", bin), min_max |
| 157 | +end |
0 commit comments