octave_packages/image-1.0.15/qtdecomp.m

   1 ## Copyright (C) 2004 Josep Mones i Teixidor
   2 ##
   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.
   7 ##
   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.
  12 ##
  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/>.
  15
  16 ## -*- texinfo -*-
  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.
  24 ##
  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.
  29 ##
  30 ## Note that blocks are not decomposed if their dimensions are not even.
  31 ##
  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).
  36 ##
  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.
  39 ##
  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.
  46 ##
  47 ## S = qtdecomp(I, threshold, mindim) decomposes an image using the
  48 ## @var{threshold} as just described, but doesn't produce blocks smaller
  49 ## than mindim.
  50 ##
  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}
  54 ## was considered.
  55 ##
  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
  62 ## otherwise.
  63 ##
  64 ## S = qtdecomp(I, fun, ...) behaves as qtdecomp(I, fun) but passes
  65 ## extra parameters to @var{fun}.
  66 ##
  67 ## @seealso{qtgetblk, qtsetblk}
  68 ## @end deftypefn
  69
  70 ## Author:  Josep Mones i Teixidor <jmones@puntbarra.com>
  71
  72 function S = qtdecomp(I, p1, varargin)
  73   if (nargin<1)
  74     usage("S=qtdecomp(I)");
  75   endif
  76
  77   if (!ismatrix(I) || size(I,1)!=size(I,2))
  78     error("qtdecomp: I should be a square matrix.");
  79   endif
  80
  81   ## current size (assumed to be square)
  82   curr_size=size(I,1);
  83
  84   ## initial mindim to a sensible value
  85   mindim=1;
  86
  87   ## sensible default maxdim value
  88   maxdim=curr_size;
  89
  90   if (nargin<2)
  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
  96     ## expensive.
  97     decision_method=0;
  98   elseif (isreal(p1))
  99     ## p1 is threshold
 100     threshold=p1;
 101     decision_method=1;
 102
 103     if(strcmp(typeinfo(I), 'uint8 matrix'))
 104       threshold*=255;
 105     elseif(strcmp(typeinfo(I), 'uint16 matrix'))
 106       threshold*=65535;
 107     endif
 108
 109     if (nargin>3)
 110       usage("S=qtdecomp(I,threshold,mindim),        \
 111           S=qtdecomp(I,threshold,[mindim maxdim])");
 112     elseif (nargin==3)
 113       dims=varargin{1};
 114       if (isvector(dims)&&length(dims)==2)
 115         mindim=dims(1);
 116         maxdim=dims(2);
 117       elseif (isreal(dims))
 118         mindim=dims;
 119       else
 120         error("qtdecomp: third parameter must be 'mindim' or '[mindim maxdim]'");
 121       endif
 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.
 124     endif
 125
 126   elseif strcmp(typeinfo(p1),"function handle") ...
 127           || strcmp(typeinfo(p1),"inline function")
 128     ## function handles seem to return true to isscalar
 129     fun=p1;
 130     decision_method=2;
 131   else
 132     error("qtdecomp: second parameter must be a integer (threshold) or a function handle (fun).");
 133   endif
 134
 135   ## initialize results matrices
 136   res=[];
 137
 138   ## bool to flag end
 139   finished=false;
 140
 141   ## array of offsets to blocks to evaluate
 142   offsets=[1,1];
 143
 144   if (maxdim<mindim)
 145     error("qtdecomp: mindim must be smaller than maxdim.");
 146   endif
 147
 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));
 151     if(initial_splits>0)
 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.");
 155       endif
 156       ## update curr_size
 157       curr_size/=divs;
 158       if(curr_size<mindim)
 159         error("qtdecomp: maxdim restriction collides with mindim restriction.");
 160       endif
 161       els=([0:divs-1]*curr_size+1).';
 162       offsets=[kron(els,ones(divs,1)), kron(ones(divs,1),els)];
 163     endif
 164   endif
 165
 166   while(!finished && rows(offsets)>0)
 167     ## check other ending conditions:
 168     ## is size is odd?
 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];
 173       finished = true;
 174     else
 175       if (decision_method<2)
 176         db=logical(ones(rows(offsets),1));
 177         for r=1:rows(offsets)
 178           o=offsets(r,:);
 179           fo=offsets(r,:)+curr_size-1;
 180
 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))))
 184               db(r)=0;
 185             endif
 186           else
 187             ## check threshold
 188             t=I(o(1):fo(1),o(2):fo(2));
 189             t=t(:);
 190             if ((max(t)-min(t))<=threshold)
 191               db(r)=0;
 192             endif
 193           endif
 194         endfor
 195       elseif(decision_method==2)
 196         ## function handle decision method
 197         ## build blocks
 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));
 202         endfor
 203
 204         db=feval(fun, b, varargin{:});
 205       else
 206         error("qtdecomp: execution shouldn't reach here. Please report this as a bug.");
 207       endif
 208
 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];
 212
 213       ## Update curr_size for next iteration
 214       curr_size/=2;
 215
 216       ## Prepare offsets for next iteration
 217       otemp=offsets(find(db),:);
 218       hs=ones(rows(otemp),1)*curr_size;
 219       zs=zeros(size(hs));
 220       offsets=[otemp;otemp+[hs,zs];otemp+[zs,hs];otemp+[hs,hs]];
 221     endif
 222   endwhile
 223
 224   S=sparse(res(:,1),res(:,2),res(:,3),size(I,1),size(I,2));
 225 endfunction
 226
 227
 228 %!demo
 229 %! full(qtdecomp(eye(8)))
 230 %! %It finds 2 big blocks of 0 and it decomposes further where 0 and 1 are mixed.
 231
 232
 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));
 236
 237 %!# Test 'equal' method
 238 %!test
 239 %! a=ones(2,2);
 240 %! b=[2,0;0,0];
 241 %! assert(full(qtdecomp(eye(4))), [a,b;b,a]);
 242
 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];
 252 %! B2=[2,0;0,0];
 253 %! B4=zeros(4); B4(1,1)=4;
 254
 255 %!test
 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);
 259
 260 %!# Test 'threshold' method
 261 %!test
 262 %! R=[ones(4,8); [ones(2),B2;B2,ones(2)],B4];
 263 %! assert(full(qtdecomp(A,1)), R);
 264
 265 %!test
 266 %! R=[[B4,[B2,B2;B2,B2]]; [[ones(2),B2;B2,B2],B4]];
 267 %! assert(full(qtdecomp(A,10)), R);
 268
 269 %!test
 270 %! R=[[B4,[B2,B2;B2,B2]]; [[B2,B2;B2,B2],B4]];
 271 %! assert(full(qtdecomp(A,10,2)), R);
 272 %!
 273 %! assert(full(qtdecomp(A,100,[2, 4])), [B4,B4;B4,B4]);
 274
 275 %!test
 276 %! f = @(A, c1 = 54, c2 = 0, c3 = 0) y = (A (1, 1, :) != ((c1+c2+c3) * ones (1, 1, size (A, 3))))(:);
 277 %!
 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)]);
 281
 282 %!test
 283 %!# no params
 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)]);
 286
 287 %!test
 288 %!# 1 param
 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)]);
 291
 292 %!test
 293 %!# 3 params
 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)]);
 296
 297
 298
 299 %
 300 % $Log$
 301 % Revision 1.5  2007/03/23 16:14:37  adb014
 302 % Update the FSF address
 303 %
 304 % Revision 1.4  2007/01/04 23:50:47  hauberg
 305 % Put seealso before end deftypefn
 306 %
 307 % Revision 1.3  2007/01/04 23:41:47  hauberg
 308 % Minor changes in help text
 309 %
 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
 312 %
 313 % Revision 1.1  2006/08/20 12:59:35  hauberg
 314 % Changed the structure to match the package system
 315 %
 316 % Revision 1.8  2006/01/03 02:09:15  pkienzle
 317 % Reorder tests so that shared variables are displayed unless relevant.
 318 %
 319 % Revision 1.7  2006/01/02 22:05:06  pkienzle
 320 % Reduce number of shared variables in tests
 321 %
 322 % Revision 1.6  2004/09/09 19:36:35  jmones
 323 % all_va_args -> varargin{:}. Now works on 2.1.58
 324 %
 325 % Revision 1.5  2004/09/08 14:07:22  pkienzle
 326 % Fix test for 'inline function'
 327 %
 328 % Revision 1.4  2004/08/11 19:52:41  jmones
 329 % qtsetblk added
 330 %
 331 % Revision 1.3  2004/08/11 00:05:21  jmones
 332 % seealso qtgetblk added to doc
 333 %
 334 % Revision 1.2  2004/08/10 00:19:42  jmones
 335 % Corrected misleading comment.
 336 %
 337 % Revision 1.1  2004/08/09 01:48:54  jmones
 338 % Added qtdecomp: quadtree decomposition
 339 %
 340 %
 341