1## Copyright (C) 2013 Carnë Draug <carandraug@octave.org> 2## 3## This program is free software; you can redistribute it and/or modify it under 4## the terms of the GNU General Public License as published by the Free Software 5## Foundation; either version 3 of the License, or (at your option) any later 6## version. 7## 8## This program is distributed in the hope that it will be useful, but WITHOUT 9## ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or 10## FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more 11## details. 12## 13## You should have received a copy of the GNU General Public License along with 14## this program; if not, see <http://www.gnu.org/licenses/>. 15 16## -*- texinfo -*- 17## @deftypefn {Function File} {} im2col (@var{A}, @var{block_size}) 18## @deftypefnx {Function File} {} im2col (@var{A}, @var{block_size}, @var{block_type}) 19## @deftypefnx {Function File} {} im2col (@var{A}, "indexed", @dots{}) 20## Rearrange blocks from matrix into columns. 21## 22## Rearranges blocks of size @var{block_size}, sampled from the matrix @var{A}, 23## into a serie of columns. This effectively transforms any image into a 24## 2 dimensional matrix, a block per column, which can then be passed to 25## other functions that perform calculations along columns. 26## 27## Both blocks and matrix @var{A} can have any number of dimensions (though 28## for sliding blocks, a block can't be larger than @var{A} in any dimension). 29## Blocks are always accessed in column-major order (like Octave arrays are 30## stored) so that a matrix can be easily reconstructed with @code{reshape} 31## and @code{col2im}. For a 2 dimensional matrix, blocks are taken first from 32## the top to the bottom, and then from the left to the right of the matrix. 33## 34## The sampling can be performed in two different ways as defined by 35## @var{block_type} (defaults to @qcode{"sliding"}): 36## 37## @table @asis 38## @item @qcode{"distinct"} 39## Each block is completely distinct from the other, with no overlapping 40## elements. The matrix @var{A} is padded as required with a value of 0 41## (or 1 for non-integer indexed images). 42## 43## @item @qcode{"sliding"} 44## A single block slides across @var{A} without any padding. 45## 46## While this can be used to perform sliding window operations such as maximum 47## and median filters, specialized functions such as @code{imdilate} and 48## @code{medfilt2} will be more efficient. 49## 50## Note that large images being arranged in large blocks can easily exceed the 51## maximum matrix size (see @code{sizemax}). For example, a matrix @var{A} of 52## size 500x500, with sliding block of size [100 100], would require a matrix 53## with 2.4108e+09 elements, i.e., the number of elements in a block, 54## @code{100*100}, times the number of blocks, @code{(500-10+1) * (500-10+1)}. 55## 56## @end table 57## 58## If @var{A} is an indexed image, the second argument should be the 59## string @qcode{"indexed"} so that any required padding is done correctly. 60## The padding value will be 0 except for indexed images of class uint8 61## and uint16. 62## 63## @seealso{blockproc, bestblk, col2im, colfilt, nlfilter, reshape} 64## @end deftypefn 65 66## Matlab behaves weird with N-dimensional images. It ignores block_size 67## elements after the first 2, and treat N-dimensional as if the extra 68## dimensions were concatenated horizontally. We are performing real 69## N-dimensional conversion of image blocks into colums. 70 71function B = im2col (A, varargin) 72 73 ## Input check 74 if (nargin > 4) 75 print_usage (); 76 endif 77 [p, block_size, padval] = im2col_check ("im2col", nargin, A, varargin{:}); 78 if (nargin > p) 79 ## we have block_type param 80 if (! ischar (varargin{p})) 81 error("im2col: BLOCK_TYPE must be a string"); 82 endif 83 block_type = varargin{p++}; 84 else 85 block_type = "sliding"; 86 endif 87 if (nargin > p) 88 print_usage (); 89 endif 90 91 ## After all the input check, start the actual im2col. The idea is to 92 ## calculate the linear indices for each of the blocks (using broadcasting 93 ## for each dimension), and then reshape it with one block per column. 94 95 switch (tolower (block_type)) 96 case "distinct" 97 ## We may need to expand the size vector to include singletons 98 size_singletons = @(x, ndim) postpad (size (x), ndim, 1); 99 100 ## Calculate needed padding 101 A_size = size_singletons (A, numel (block_size)); 102 sp = mod (-A_size, block_size); 103 if (any (sp)) 104 A = padarray (A, sp, padval, "post"); 105 endif 106 A_size = size_singletons (A, numel (block_size)); 107 108 ## Get linear indixes for the first block 109 [ind, stride] = get_1st_ind (A_size, block_size); 110 111 ## Get linear indices for all of the blocks 112 blocks = A_size ./ block_size; 113 step = stride .* block_size; 114 limit = step .* (blocks -1); 115 for dim = 1:numel (A_size) 116 ind = ind(:) .+ (0:step(dim):limit(dim)); 117 endfor 118 n_blocks = prod (blocks); 119 120 case "sliding" 121 if (numel (block_size) > ndims (A)) 122 error ("im2col: BLOCK_SIZE can't have more elements than the dimensions of A"); 123 endif 124 125 ## Get linear indixes for the first block 126 [ind, stride] = get_1st_ind (size (A), block_size); 127 128 ## Get linear indices for all of the blocks 129 slides = size (A) - block_size; 130 limit = stride .* slides; 131 for dim = 1:ndims (A) 132 ## We need to use bsxfun here because of 133 ## https://savannah.gnu.org/bugs/?47085 134 ind = bsxfun (@plus, ind(:), (0:stride(dim):limit(dim))); 135 endfor 136 n_blocks = prod (slides +1); 137 138 otherwise 139 error ("im2col: invalid BLOCK_TYPE `%s'.", block_type); 140 endswitch 141 142 B = reshape (A(ind(:)), prod (block_size), n_blocks); 143endfunction 144 145## Get linear indices and for the first block, and stride size per dimension 146function [ind, stride] = get_1st_ind (A_size, block_size) 147 stride = [1 cumprod(A_size(1:end-1))]; 148 limit = (block_size -1) .* stride; 149 ind = 1; 150 for dim = 1:numel (A_size) 151 ind = ind(:) .+ (0:stride(dim):limit(dim)); 152 endfor 153endfunction 154 155%!demo 156%! ## Divide A using distinct blocks and then reverse the operation 157%! A = [ 1:10 158%! 11:20 159%! 21:30 160%! 31:40]; 161%! B = im2col (A, [2 5], "distinct") 162%! C = col2im (B, [2 5], [4 10], "distinct") 163 164## test default block type 165%!test 166%! a = rand (10); 167%! assert (im2col (a, [5 5]), im2col (a, [5 5], "sliding")) 168 169## indexed makes no difference when sliding 170%!test 171%! a = rand (10); 172%! assert (im2col (a, [5 5]), im2col (a, "indexed", [5 5])) 173 174%!error <BLOCK_TYPE> im2col (rand (20), [2 5], 10) 175%!error <BLOCK_TYPE> im2col (rand (20), [2 5], "wrong_block_type") 176%!error im2col (rand (10), [5 5], "sliding", 5) 177%!error im2col (rand (10), "indexed", [5 5], "sliding", 5) 178 179%!shared B, A, Bs, As, Ap, Bp0, Bp1, Bp0_3s 180%! v = [1:10]'; 181%! r = reshape (v, 2, 5); 182%! B = [v v+20 v+40 v+10 v+30 v+50]; 183%! A = [r r+10; r+20 r+30; r+40 r+50]; 184%! As = [ 1 2 3 4 5 185%! 6 7 8 9 10 186%! 11 12 13 14 15]; 187%! b1 = As(1:2, 1:4)(:); 188%! b2 = As(2:3, 1:4)(:); 189%! b3 = As(1:2, 2:5)(:); 190%! b4 = As(2:3, 2:5)(:); 191%! Bs = [b1, b2, b3, b4]; 192%! Ap = A(:, 1:9); 193%! Bp1 = Bp0 = B; 194%! Bp0(9:10, 4:6) = 0; 195%! Bp1(9:10, 4:6) = 1; 196%! Bp0_3s = Bp0; 197%! Bp0_3s(11:30, :) = 0; 198 199## test distinct block type 200%!assert (im2col (A, [2 5], "distinct"), B); 201 202## padding for distinct 203%!assert (im2col (Ap, [2 5], "distinct"), Bp0); 204%!assert (im2col (Ap, [2 5 3], "distinct"), Bp0_3s); 205%!assert (im2col (Ap, "indexed", [2 5], "distinct"), Bp1); 206%!assert (im2col (uint8 (Ap), "indexed", [2 5], "distinct"), uint8 (Bp0)); 207%!assert (im2col (uint16 (Ap), "indexed", [2 5], "distinct"), uint16 (Bp0)); 208%!assert (im2col (int16 (Ap), "indexed", [2 5], "distinct"), int16 (Bp1)); 209%!assert (im2col (uint32 (Ap), "indexed", [2 5], "distinct"), uint32 (Bp1)); 210 211## Always return correct class 212%!assert (im2col (uint8 (A), [2 5], "distinct"), uint8 (B)); 213%!assert (im2col (single (A), [2 5], "distinct"), single (B)); 214%!assert (im2col (logical (A), [2 5], "distinct"), logical (B)); 215%!assert (im2col (uint8 (As), [2 4], "sliding"), uint8 (Bs)); 216%!assert (im2col (single (As), [2 4], "sliding"), single (Bs)); 217%!assert (im2col (logical (As), [2 4], "sliding"), logical (Bs)); 218 219## test sliding block type 220%!assert (im2col (As, [2 4], "sliding"), Bs); 221%!assert (im2col (As, [3 5], "sliding"), As(:)); 222 223## Test N-dimensional 224%!test 225%! A = randi (9, 10, 9, 5); 226%!assert (convn (A, ones (3, 3, 3), "valid"), 227%! reshape (sum (im2col (A, [3 3 3])), [8 7 3])); 228%! 229%! A = randi (9, 10, 9, 5, 7); 230%!assert (convn (A, ones (3, 3, 3), "valid"), 231%! reshape (sum (im2col (A, [3 3 3])), [8 7 3 7])); 232%!assert (convn (A, ones (3, 4, 3), "valid"), 233%! reshape (sum (im2col (A, [3 4 3])), [8 6 3 7])); 234%!assert (convn (A, ones (3, 5, 3, 2), "valid"), 235%! reshape (sum (im2col (A, [3 5 3 2])), [8 5 3 6])); 236 237## Corner case for Matlab compatibility -- bug #46774 238%!assert (im2col (1:8, [2 1]), zeros (2, 0)) 239