-
Notifications
You must be signed in to change notification settings - Fork 110
/
prcclip.m
95 lines (81 loc) · 2.13 KB
/
prcclip.m
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
function [lims,A] = prcclip(A,prc,symmetric)
%PRCCLIP percentile clipping
%
% Syntax
%
% lims = prcclip(A,prc)
% lims = prcclip(A,prc,symmetric)
% [lims,Ac] = prcclip(A,prc)
%
% Description
%
% PRCCLIP determines the values lims that clip the values in A to range
% between the nth percentile and 100-nth percentile. prc can be a
% scalar or two element vector. If prc is a scalar, then lims are
% determined by the prc'th and 100-prc'th percentile. If prc is a
% vector, it contains the lower and upper percentile (e.g. [2 98]).
%
% Input arguments
%
% A GRIDobj
% prc percentile (scalar or two element vector)
% symmetric if true, lims will be symmetric around zero. prcclip will
% determine the percentiles and then will adjust the limits
% so that they range between
% [-max(abs(lims)) max(abs(lims))]
%
% Output arguments
%
% lims limits
% Ac percentile clipped GRIDobj
%
% Example
%
% DEM = GRIDobj('srtm_bigtujunga30m_utm11.tif');
% C = curvature(DEM);
% lims = prcclip(C,2,true);
% clr = ttscm('vik');
% imageschs(DEM,C,'colormap',clr,'caxis',lims);
%
%
% See also: GRIDobj/minmaxnorm
%
% Author: Wolfgang Schwanghart (w.schwanghart[at]geo.uni-potsdam.de)
% Date: 14. June, 2018
if nargin == 1
prc = 2;
symmetric = false;
elseif nargin == 2
validate
symmetric = false;
end
qclip = prc/100;
if numel(qclip) == 1
qclip = [qclip 1-qclip];
end
qclip = sort(qclip);
I = ~isnan(A.Z);
[n,edges] = histcounts(A.Z(I(:)),'Normalization','cdf');
lval = edges(find(n>=qclip(1),1,'first'));
uval = edges(find(n<(qclip(2)),1,'last'));
if lval == uval
warning('TopoToolbox:imageschs','percent clip returns flat matrix');
else
if nargout == 2
A.Z(I) = max(A.Z(I),lval);
A.Z(I) = min(A.Z(I),uval);
end
end
lims = [lval,uval];
if symmetric
lims = max(abs(lims));
lims = [-lims lims];
end
% if nargout == 2
% if flatmatrix
% A(:,:) = lval;
% else
% A = max(A,lims(1));
% A = min(A,lims(2));
% end
% end