1## Copyright (C) 2004 Josep Monés i Teixidor <jmones@puntbarra.com> 2## Copyright (C) 2015 Carnë Draug <carandraug@octave.org> 3## 4## This program is free software; you can redistribute it and/or modify 5## it under the terms of the GNU General Public License as published by 6## the Free Software Foundation; either version 3 of the License, or 7## (at your option) any later version. 8## 9## This program is distributed in the hope that it will be useful, 10## but WITHOUT ANY WARRANTY; without even the implied warranty of 11## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 12## GNU General Public License for more details. 13## 14## You should have received a copy of the GNU General Public License 15## along with this program; if not, see <http://www.gnu.org/licenses/>. 16 17## -*- texinfo -*- 18## @deftypefn {Function File} {} stretchlim (@var{I}) 19## @deftypefnx {Function File} {} stretchlim (@var{RGB}) 20## @deftypefnx {Function File} {} stretchlim (@dots{}, @var{tol}) 21## Find limits to contrast stretch an image. 22## 23## Returns a 2 element column vector, @code{[@var{low}; @var{high}]}, 24## with the pair of intensities to contrast stretch @var{I} which 25## saturates at most @var{tol} of the image. The output of this 26## function matches the input expected by @code{imadjust}. 27## 28## The input argument @var{tol}, controls the fraction of the image to be 29## saturated and defaults to @code{[0.01 0.99]}, i.e., a 1% saturation 30## on both sides of the image histogram. It can be specified in two ways: 31## 32## @itemize 33## @item a two element vector with lower and higher fraction of the 34## the image to be saturated. These values must be in the range 35## [0 1], the display range of images of floating point class. 36## 37## @item a scalar value with the fraction of image to be saturated on 38## each side; e.g., a @var{tol} with a value of @code{0.05} is equivalent 39## to @code{[0.05 0.95]}. This value must be in the range @code{[0 0.5]}. 40## 41## @end itemize 42## 43## A common use case wanting to maximize contrast without causing any 44## saturation. In such case @var{tol} should be 0 (zero). It is the 45## equivalent to @code{[min(@var{I}(:)); max(@var{I}(:))]} for a single 46## plane. 47## 48## The return values are of class double and in the range [0 1] regardless 49## of the input image class. These values are scaled in the image class 50## range (see @code{im2double}). 51## 52## If the input is a RGB image, i.e., the 3rd dimension has length 3, then 53## it returns a @code{[2 3]} matrix with separate limits for each colour. 54## It will actually do this for each plane, so an input of size, 55## @code{[M N P]} will return a @code{[2 P]} matrix. 56## 57## Higher dimensions of @var{I} are also valid. As in the 3 dimensional 58## case, the limits for each 2 dimensional plane are returned in a 2 row 59## vector for each 2 dimensional plane, one dimension below. That is, an 60## input of size @code{[M N P Q R]} will return an array of size 61## @code{[2 P Q R]}. 62## 63## Note the detail that @var{tol} is the maximum fraction of saturation. 64## It is rare that there is a value for that exact saturation. In such 65## case, @code{stretchlim} will always round down and saturate less. 66## @var{tol} is the saturation limit. For example, if @var{tol} is 67## @code{0.10}, but there are only values that will lead to 5 or 11% 68## saturation, it will return the value for a 5% saturation. 69## 70## @seealso{brighten, contrast, histeq, imadjust} 71## @end deftypefn 72 73function low_high = stretchlim (img, tol = [0.01 0.99]) 74 75 if (nargin () <1 || nargin () > 2) 76 print_usage (); 77 endif 78 79 if (! isimage (img)) 80 error("stretchlim: I or RGB must be an image"); 81 endif 82 83 ## Handle tol 84 if (nargin > 1) 85 if (! isnumeric (tol)) 86 error ("stretchlim: TOL must be numeric"); 87 endif 88 89 if (isscalar (tol)) 90 if (min (tol) < 0 || max (tol) > 0.5) 91 error ("stretchlim: TOL must be in the [0 1] range"); 92 endif 93 tol = [tol (1-tol)]; 94 95 elseif (isvector (tol)) 96 if (numel (tol) != 2) 97 error ("stretchlim: TOL must be a 2 element vector"); 98 endif 99 endif 100 endif 101 102 ## tol is about the percentage of values that will be saturated. 103 ## So instead of percentages, we convert to the actual number of 104 ## pixels that need to be saturated. After sorting the values in 105 ## the image, that number of pixels simply becomes the index for 106 ## the limits. 107 ## 108 ## Note that the actual intensity value that we set the limits to, 109 ## is not saturated. Only the values below or above the lower and 110 ## higher limits it will be considered saturated. 111 ## 112 ## And since most images will have repeated values in the pixels, 113 ## chances are that there's not a limit that would cause only the 114 ## exact percentage of pixels to be saturated. In such cases, we 115 ## must prefer a limit that would saturate less pixels than the 116 ## requested, rather than the opposite. 117 ## 118 ## We want to compute this for each plane, so we reshape the image 119 ## in order to have each plane into a single column, while respecting 120 ## any other dimensions beyond the 3rd. 121 122 sz = postpad (size (img), max (3, ndims (img)), 1); 123 plane_length = sz(1) * sz(2); 124 125 img = reshape (img, [plane_length sz(3:end)]); 126 127 lo_idx = floor (tol(1) * plane_length) + 1; 128 hi_idx = ceil (tol(2) * plane_length); 129 130 if (lo_idx == 1 && hi_idx == plane_length) 131 ## special case, equivalent to tol [0 1], even if tol was not 132 ## actually [0 1] but the image size would effectively make it. 133 low_high = [min(img, [], 1); max(img, [], 1)]; 134 else 135 strides = reshape ((0:plane_length:(numel(img)-1)), [1 sz(3:end)]); 136 lo_hi_idx = [lo_idx; hi_idx] .+ strides; 137 sorted = sort (img, 1); 138 low_high = sorted(lo_hi_idx); 139 endif 140 141 low_high = im2double (low_high); 142 143 ## This will only happen if the input was floating point and with 144 ## values already outside the [0 1] range (common for users to simply 145 ## convert an uint8 image to double and forget that display is [0 1] 146 ## only). 147 low_high(low_high < 0) = 0; 148 low_high(low_high > 1) = 1; 149 150 ## If a plane has a single value, then both limits will have the same 151 ## value. In such case, we must return [0; 1] 152 equal_cols = low_high(1,:) == low_high(2,:); 153 low_high(:, equal_cols) = repmat ([0; 1], 1, nnz (equal_cols)); 154endfunction 155 156%!error (stretchlim ()); 157%!error (stretchlim ("bad parameter")); 158#%!error (stretchlim (zeros (10, 10, 3, 2))); 159%!error (stretchlim (zeros (10, 10), "bad parameter")); 160%!error (stretchlim (zeros (10, 10), 0.01, 2)); 161 162## default parameters 163%!assert (stretchlim (0.01:.01:1), [0.02; 0.99]) 164%!assert (stretchlim (0.01:.01:1), stretchlim (0.01:.01:1, [0.01 0.99])) 165 166## use scalar for tol 167%!assert (stretchlim (0.01:.01:1, 0.15), stretchlim (0.01:.01:1, [0.15 0.85])) 168 169## this is different than Matlab but it looks like it's a Matlab bug 170## (Matlab returns [0.018997482261387; 0.951003280689708]) 171## We actually have differences from Matlab which sometimes returns 172## values that are not present in the image. 173%!assert (stretchlim (0.01:.01:1, [0.01,0.95]), [0.02; 0.95], eps) 174 175## corner case of zero tolerance 176%!assert (stretchlim (0.01:.01:1, 0), [0.01; 1]) 177 178%!test 179%! im = rand (5); 180%! assert (stretchlim (im, 0), [min(im(:)); max(im(:))]) 181 182%!test 183%! im = rand (5, 5, 3); 184%! assert (stretchlim (im, 0), 185%! [min(im(:,:,1)(:)) min(im(:,:,2)(:)) min(im(:,:,3)(:)); 186%! max(im(:,:,1)(:)) max(im(:,:,2)(:)) max(im(:,:,3)(:))]) 187 188 189## corner case where tol is not zero but the image is so small that 190## it might as well be. 191%!test 192%! im = rand (5); 193%! assert (stretchlim (im, 0.03), [min(im(:)); max(im(:))]) 194%! assert (stretchlim (im, 0.0399), [min(im(:)); max(im(:))]) 195 196 197## Test with non double data-types 198%!assert (stretchlim (uint8 (1:100)), im2double (uint8 ([2; 99]))) 199%!assert (stretchlim (uint8 (1:100), .25), im2double (uint8 ([26; 75]))) 200%!assert (stretchlim (uint16 (1:1000)), im2double (uint16 ([11; 990]))) 201 202%!assert (stretchlim (int16 (-100:100)), im2double (int16 ([-98; 98]))) 203%!assert (stretchlim (single (0.01:.01:1)), 204%! double (single (0.01:.01:1)([2; 99])).') 205 206 207## non uniform histogram tests 208%!assert (stretchlim (uint8 ([1 repmat(2, [1, 90]) 92:100]), 0.05), 209%! im2double (uint8 ([2; 95]))) 210%!assert (stretchlim (uint8 ([1 repmat(2, [1 4]) 6:100]), 0.05), 211%! im2double (uint8 ([6; 95]))) 212 213## test limit rounding (actually, lack of rounding, we always round down) 214## Note that this tests were different in the image package before v2.6. 215## Back then we performed rounding of the fraction that was saturated. 216%!assert (stretchlim (uint8 ([1 repmat(2, [1 5]) 7:100]), 0.05), 217%! im2double (uint8 ([2; 95]))) 218%!assert (stretchlim (uint8 ([1 repmat(2, [1 6]) 8:100]), 0.05), 219%! im2double (uint8 ([2; 95]))) 220%!assert (stretchlim (uint8 ([1 repmat(2, [1 7]) 9:100]), 0.05), 221%! im2double (uint8 ([2; 95]))) 222%!assert (stretchlim (uint8 ([1 repmat(2, [1 8]) 10:100]), 0.05), 223%! im2double (uint8 ([2; 95]))) 224 225%!assert (stretchlim (uint8 ([1 repmat(2, [1 5]) repmat(3, [1 5]) 9:100]), 0.04), 226%! im2double (uint8 ([2; 96]))) 227%!assert (stretchlim (uint8 ([1 repmat(2, [1 5]) repmat(3, [1 5]) 9:100]), 0.05), 228%! im2double (uint8 ([2; 95]))) 229%!assert (stretchlim (uint8 ([1 repmat(2, [1 5]) repmat(3, [1 5]) 9:100]), 0.06), 230%! im2double (uint8 ([3; 94]))) 231%!assert (stretchlim (uint8 ([1 repmat(2, [1 5]) repmat(3, [1 5]) 9:100]), 0.07), 232%! im2double (uint8 ([3; 93]))) 233%!assert (stretchlim (uint8 ([1 repmat(2, [1 5]) repmat(3, [1 5]) 9:100]), 0.08), 234%! im2double (uint8 ([3; 92]))) 235 236## test RGB 237%!test 238%! RGB = zeros (100, 1, 3, "uint16"); 239%! RGB(:,:,1) = [1:1:100]; 240%! RGB(:,:,2) = [2:2:200]; 241%! RGB(:,:,3) = [4:4:400]; 242%! assert (stretchlim (RGB) , im2double (uint16 ([2 4 8; 99 198 396]))) 243 244## test other 3D lengths 245%!test 246%! im6c = zeros (100, 1, 6, "uint16"); 247%! im6c(:,:,1) = [1:1:100]; 248%! im6c(:,:,2) = [2:2:200]; 249%! im6c(:,:,3) = [4:4:400]; 250%! im6c(:,:,4) = [8:8:800]; 251%! im6c(:,:,5) = [16:16:1600]; 252%! im6c(:,:,6) = [32:32:3200]; 253%! assert (stretchlim (im6c) , 254%! im2double (uint16 ([2 4 8 16 32 64; 99 198 396 792 1584 3168]))) 255 256%!test 257%! im = [0 0 .1 .1 .1 .1 .2 .2 .2 .4 .4 .6 .6 .7 .7 .9 .9 .9 1 1]; 258%! 259%! assert (stretchlim (im), [0; 1]) 260%! 261%! ## Consider the returned lower limit in this test. A lower limit 262%! ## of 0.1 will saturate two elements (10%), while 0.2 will saturate 263%! ## 6 elements (30%). Both have the same distance to 20% but returning 264%! ## 0.1 is Matlab compatible. 265%! ## Now looking at the higher limit. A limit of .9 will saturate 266%! ## 2 elements (10%), while a limit of 0.7 will saturate 5 elements (25%). 267%! ## However, for Matlab compatibility we must return .9 even though 268%! ## 25% would be closer to 20%. 269%! ## Basically, it's not just rounded. 270%! assert (stretchlim (im, .2), [0.1; 0.9]) 271%! 272%! assert (stretchlim (im, .15), [0.1; 0.9]) 273%! assert (stretchlim (im, .1), [0.1; 0.9]) 274%! assert (stretchlim (im, .25), [0.1; 0.7]) 275%! 276%! ## Reorder the vector of values (real images don't have the values 277%! ## already sorted), just to be sure it all works. 278%! im([6 3 16 11 7 17 14 8 5 19 15 1 2 4 18 13 9 20 10 12]) = im; 279%! assert (stretchlim (im, .2), [0.1; 0.9]) 280%! assert (stretchlim (im, .15), [0.1; 0.9]) 281%! assert (stretchlim (im, .1), [0.1; 0.9]) 282%! assert (stretchlim (im, .25), [0.1; 0.7]) 283 284## odd length images to test rounding of saturated fraction. With a 1% 285## fraction to be saturated and 991 elements, that's 9.91 pixels. Since 286## TOL is the limit, we must saturate the top and bottom 9 pixels (not 10). 287%!assert (stretchlim (0.01:.001:1), [0.019; 0.991], eps) 288%!assert (stretchlim (0.01:.001:1, [0.01,0.95]), [0.019; 0.951], eps) 289%!assert (stretchlim (0.01:.001:1, 0), [0.01; 1]) 290%!assert (stretchlim (single (0.01:.001:1)), 291%! double (single (0.01:.001:1)([10; 982])).') 292 293## Test ignoring floating point values outside [0 1] and how they do 294## count for the fraction of saturated values, but are simply clipped 295## to the [0 1] range 296%!test 297%! assert (stretchlim ([(.05:.05:1) (2:4)], 0.2), [0.25; 0.95], eps) 298%! assert (stretchlim ([(.05:.05:1) (2:5)], 0.2), [0.25; 1]) 299%! assert (stretchlim ([(.05:.05:1) (2:6)], 0.2), [0.3; 1]) 300%! assert (stretchlim ([(.05:.05:1) (2:7)], 0.2), [0.3; 1]) 301 302%!test 303%! assert (stretchlim ([(-6:0) (.05:.05:1)], 0.2), [0; 0.75], eps) 304%! assert (stretchlim ([(-5:0) (.05:.05:1)], 0.2), [0; 0.75], eps) 305 306## Test N dimensions 307%!test 308%! im = rand (4, 4, 2, 3, 2); 309%! rv = zeros (2, 2, 3, 2); 310%! for p = 1:2 311%! for q = 1:3 312%! for r = 1:2 313%! rv(:,p,q,r) = stretchlim (im(:,:,p,q,r), 0.25); 314%! endfor 315%! endfor 316%! endfor 317%! assert (stretchlim (im, 0.25), rv) 318 319## The same but important because of our special code path for tol == 0 320%!test 321%! im = rand (4, 4, 2, 3, 2); 322%! rv = zeros (2, 2, 3, 2); 323%! for p = 1:2 324%! for q = 1:3 325%! for r = 1:2 326%! rv(:,p,q,r) = stretchlim (im(:,:,p,q,r), 0); 327%! endfor 328%! endfor 329%! endfor 330%! assert (stretchlim (im, 0), rv) 331 332## Corner case there's a single value in the whole image 333%!assert (stretchlim (zeros (5)), [0; 1]) 334%!assert (stretchlim (ones (5)), [0; 1]) 335%!assert (stretchlim (.6 * ones (5)), [0; 1]) 336%!assert (stretchlim (zeros (3, 3, 3, 3)), repmat ([0; 1], [1 3 3])) 337%!assert (stretchlim ([0 .5 .5 .5 .5 1], .2), [0; 1]) 338 339## Cornier case when some of the planes have a single unique value 340%!test 341%! im = repmat ((magic (5) -1) / 24, [1 1 3 3]); 342%! im(:,:,1,1) = 0; 343%! im(:,:,2,2) = .5; 344%! im(:,:,3,3) = 1; 345%! lims = stretchlim (im, 0.2); 346%! assert (size (lims), [2 3 3]) 347%! assert (lims(:, [2 3 4 6 7 8]), 348%! repmat ([(1/24)*round(24*.2); 1-((1/24)*round(24*.2))], [1 6]), eps) 349%! assert (lims(:, [1 5 9]), repmat ([0; 1], [1 3])) 350