1 ## Copyright (C) 2004 Josep Mones i Teixidor
3 ## This program is free software; you can redistribute it and/or modify
4 ## it under the terms of the GNU General Public License as published by
5 ## the Free Software Foundation; either version 2 of the License, or
6 ## (at your option) any later version.
8 ## This program is distributed in the hope that it will be useful,
9 ## but WITHOUT ANY WARRANTY; without even the implied warranty of
10 ## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
11 ## GNU General Public License for more details.
13 ## You should have received a copy of the GNU General Public License
14 ## along with this program; If not, see <http://www.gnu.org/licenses/>.
17 ## @deftypefn {Function File} {@var{S} = } qtdecomp (@var{I})
18 ## @deftypefnx {Function File} {@var{S} = } qtdecomp (@var{I},@var{threshold})
19 ## @deftypefnx {Function File} {@var{S} = } qtdecomp (@var{I},@var{threshold},@var{mindim})
20 ## @deftypefnx {Function File} {@var{S} = } qtdecomp (@var{I},@var{threshold},@var{[mindim maxdim]})
21 ## @deftypefnx {Function File} {@var{S} = } qtdecomp (@var{I},@var{fun})
22 ## @deftypefnx {Function File} {@var{S} = } qtdecomp (@var{I},@var{fun},@var{P1},@var{P2},...)
23 ## Performs quadtree decomposition.
25 ## qtdecomp decomposes a square image @var{I} into four equal-sized
26 ## blocks. Then it performs some kind of test on each block to decide if
27 ## it should decompose them further. This process is repeated
28 ## iteratively until there's no block left to be decomposed.
30 ## Note that blocks are not decomposed if their dimensions are not even.
32 ## The output is a sparse matrix whose non-zero elements determine the
33 ## position of the block (the element is at top-left position in the
34 ## block) and size of each block (the value of the element determines
35 ## length of a side of the square-shaped block).
37 ## S = qtdecomp(I) decomposes an intensity image @var{I} as described
38 ## above. By default it doesn't split a block if all elements are equal.
40 ## S = qtdecomp(I, threshold) decomposes an image as decribed, but only
41 ## splits a block if the maximum value in the block minus the minimum
42 ## value is greater than @var{threshold}, which is a value between 0 and
43 ## 1. If @var{I} is of class uint8, @var{threshold} is multiplied by 255
44 ## before use. Also, if@var{I} is of class uint16, @var{threshold} is
45 ## multiplied by 65535.
47 ## S = qtdecomp(I, threshold, mindim) decomposes an image using the
48 ## @var{threshold} as just described, but doesn't produce blocks smaller
51 ## S = qtdecomp(I, threshold, [mindim maxdim]) decomposes an image as
52 ## described, but produces blocks that can't be bigger than maxdim. It
53 ## decomposes to maxdim even if it isn't needed if only @var{threshold}
56 ## S = qtdecomp(I, fun) decomposes an image @var{I} and uses function
57 ## @var{fun} to decide if a block should be splitted or not. @var{fun}
58 ## is called with a m-by-m-by-k array of m-by-m blocks to be
59 ## considered, and should return a vector of size k, whose elements
60 ## represent each block in the stacked array. @var{fun} sets the
61 ## corresponding value to 1 if the block should be split, and 0
64 ## S = qtdecomp(I, fun, ...) behaves as qtdecomp(I, fun) but passes
65 ## extra parameters to @var{fun}.
67 ## @seealso{qtgetblk, qtsetblk}
70 ## Author: Josep Mones i Teixidor <jmones@puntbarra.com>
72 function S = qtdecomp(I, p1, varargin)
74 usage("S=qtdecomp(I)");
77 if (!ismatrix(I) || size(I,1)!=size(I,2))
78 error("qtdecomp: I should be a square matrix.");
81 ## current size (assumed to be square)
84 ## initial mindim to a sensible value
87 ## sensible default maxdim value
91 ## Initialize decision method variable
92 ## We could have implemented threshold as a function and use an
93 ## uniform interface (function handle) to decide whether to split or
94 ## not blocks. We have decided not to do so because block
95 ## rearrangement that is needed as a parameter to functions is
103 if(strcmp(typeinfo(I), 'uint8 matrix'))
105 elseif(strcmp(typeinfo(I), 'uint16 matrix'))
110 usage("S=qtdecomp(I,threshold,mindim), \
111 S=qtdecomp(I,threshold,[mindim maxdim])");
114 if (isvector(dims)&&length(dims)==2)
117 elseif (isreal(dims))
120 error("qtdecomp: third parameter must be 'mindim' or '[mindim maxdim]'");
122 ## we won't check if mindim or maxdim are powers of 2. It's too
123 ## restrictive and don't need it at all.
126 elseif strcmp(typeinfo(p1),"function handle") ...
127 || strcmp(typeinfo(p1),"inline function")
128 ## function handles seem to return true to isscalar
132 error("qtdecomp: second parameter must be a integer (threshold) or a function handle (fun).");
135 ## initialize results matrices
141 ## array of offsets to blocks to evaluate
145 error("qtdecomp: mindim must be smaller than maxdim.");
148 ## See if we have to split a minimum regarless other considerations.
149 if (maxdim<curr_size)
150 initial_splits=ceil(log2(curr_size/maxdim));
152 divs=2^initial_splits;
153 if (rem(curr_size,divs)!=0)
154 error("qtdecomp: Can't decompose I enough times to fulfill maxdim requirement.");
159 error("qtdecomp: maxdim restriction collides with mindim restriction.");
161 els=([0:divs-1]*curr_size+1).';
162 offsets=[kron(els,ones(divs,1)), kron(ones(divs,1),els)];
166 while(!finished && rows(offsets)>0)
167 ## check other ending conditions:
169 ## is splitted size < than mindim?
170 if ((rem(curr_size,2)!=0)||((curr_size/2)<mindim))
171 ## can't continue, lets add current evaluation blocks to results
172 res=[res; offsets, ones(size(offsets,1),1)*curr_size];
175 if (decision_method<2)
176 db=logical(ones(rows(offsets),1));
177 for r=1:rows(offsets)
179 fo=offsets(r,:)+curr_size-1;
181 if(decision_method==0)
182 ## is everything equal?
183 if (all(I(o(1),o(2))==I(o(1):fo(1),o(2):fo(2))))
188 t=I(o(1):fo(1),o(2):fo(2));
190 if ((max(t)-min(t))<=threshold)
195 elseif(decision_method==2)
196 ## function handle decision method
198 b=zeros(curr_size,curr_size,rows(offsets));
199 rbc=offsets(:,1:2)+curr_size-1;
200 for r=1:rows(offsets)
201 b(:,:,r)=I(offsets(r,1):rbc(r,1),offsets(r,2):rbc(r,2));
204 db=feval(fun, b, varargin{:});
206 error("qtdecomp: execution shouldn't reach here. Please report this as a bug.");
209 ## Add blocks that won't divide to results
210 nd=offsets(find(!db),:);
211 res=[res; nd, ones(size(nd,1),1)*curr_size];
213 ## Update curr_size for next iteration
216 ## Prepare offsets for next iteration
217 otemp=offsets(find(db),:);
218 hs=ones(rows(otemp),1)*curr_size;
220 offsets=[otemp;otemp+[hs,zs];otemp+[zs,hs];otemp+[hs,hs]];
224 S=sparse(res(:,1),res(:,2),res(:,3),size(I,1),size(I,2));
229 %! full(qtdecomp(eye(8)))
230 %! %It finds 2 big blocks of 0 and it decomposes further where 0 and 1 are mixed.
233 %!# Test if odd-sized limits split
234 %!assert(full(qtdecomp(eye(5))), reshape([5,zeros(1,24)],5,5));
235 %!assert(full(qtdecomp(eye(6))), repmat(reshape([3,zeros(1,8)],3,3),2,2));
237 %!# Test 'equal' method
241 %! assert(full(qtdecomp(eye(4))), [a,b;b,a]);
243 %!shared A, B2, B4, f
244 %! A=[ 1, 4, 2, 5,54,55,61,62;
245 %! 3, 6, 3, 1,58,53,67,65;
246 %! 3, 6, 3, 1,58,53,67,65;
247 %! 3, 6, 3, 1,58,53,67,65;
248 %! 23,42,42,42,99,99,99,99;
249 %! 27,42,42,42,99,99,99,99;
250 %! 23,22,26,25,99,99,99,99;
251 %! 22,22,24,22,99,99,99,99];
253 %! B4=zeros(4); B4(1,1)=4;
256 %! R=[ones(4,8); [ones(2),B2;ones(2,4)], B4];
257 %! assert(full(qtdecomp(A)), R);
258 %! assert(full(qtdecomp(A,0)), R);
260 %!# Test 'threshold' method
262 %! R=[ones(4,8); [ones(2),B2;B2,ones(2)],B4];
263 %! assert(full(qtdecomp(A,1)), R);
266 %! R=[[B4,[B2,B2;B2,B2]]; [[ones(2),B2;B2,B2],B4]];
267 %! assert(full(qtdecomp(A,10)), R);
270 %! R=[[B4,[B2,B2;B2,B2]]; [[B2,B2;B2,B2],B4]];
271 %! assert(full(qtdecomp(A,10,2)), R);
273 %! assert(full(qtdecomp(A,100,[2, 4])), [B4,B4;B4,B4]);
276 %! f = @(A, c1 = 54, c2 = 0, c3 = 0) y = (A (1, 1, :) != ((c1+c2+c3) * ones (1, 1, size (A, 3))))(:);
278 %! assert(full(qtdecomp(A,f)),[ones(4),B4;ones(4,8)]);
279 %! assert(full(qtdecomp(A,f,54)),[ones(4),B4;ones(4,8)]);
280 %! assert(full(qtdecomp(A,f,4,40,10)),[ones(4),B4;ones(4,8)]);
284 %! first_eq=inline("(A(1,1,:)!=(54*ones(1,1,size(A,3))))(:)","A");
285 %! assert(full(qtdecomp(A,first_eq)),[ones(4),B4;ones(4,8)]);
289 %! first_eq=inline("(A(1,1,:)!=(c*ones(1,1,size(A,3))))(:)","A","c");
290 %! assert(full(qtdecomp(A,first_eq,54)),[ones(4),B4;ones(4,8)]);
294 %! first_eq=inline("(A(1,1,:)!=((c1+c2+c3)*ones(1,1,size(A,3))))(:)","A","c1","c2","c3");
295 %! assert(full(qtdecomp(A,first_eq,4,40,10)),[ones(4),B4;ones(4,8)]);
301 % Revision 1.5 2007/03/23 16:14:37 adb014
302 % Update the FSF address
304 % Revision 1.4 2007/01/04 23:50:47 hauberg
305 % Put seealso before end deftypefn
307 % Revision 1.3 2007/01/04 23:41:47 hauberg
308 % Minor changes in help text
310 % Revision 1.2 2006/10/09 19:58:09 adb014
311 % Remove dependency on miscellaneous transpose function. Simplify tests and add function handle tests
313 % Revision 1.1 2006/08/20 12:59:35 hauberg
314 % Changed the structure to match the package system
316 % Revision 1.8 2006/01/03 02:09:15 pkienzle
317 % Reorder tests so that shared variables are displayed unless relevant.
319 % Revision 1.7 2006/01/02 22:05:06 pkienzle
320 % Reduce number of shared variables in tests
322 % Revision 1.6 2004/09/09 19:36:35 jmones
323 % all_va_args -> varargin{:}. Now works on 2.1.58
325 % Revision 1.5 2004/09/08 14:07:22 pkienzle
326 % Fix test for 'inline function'
328 % Revision 1.4 2004/08/11 19:52:41 jmones
331 % Revision 1.3 2004/08/11 00:05:21 jmones
332 % seealso qtgetblk added to doc
334 % Revision 1.2 2004/08/10 00:19:42 jmones
335 % Corrected misleading comment.
337 % Revision 1.1 2004/08/09 01:48:54 jmones
338 % Added qtdecomp: quadtree decomposition