--- a/INDEX
+++ b/INDEX
@@ -18,11 +18,12 @@
  impad
  padarray 
  rotate_scale
-Statistics
+Analysis and Statistics
  corr2 
  imhist 
  mean2 
  std2 
+ qtdecomp
 Filtering
  colfilt
  histeq 

new file mode 100644
--- /dev/null
+++ b/qtdecomp.m
@@ -0,0 +1,307 @@
+## Copyright (C) 2004 Josep Mones i Teixidor
+##
+## This program is free software; you can redistribute it and/or modify
+## it under the terms of the GNU General Public License as published by
+## the Free Software Foundation; either version 2 of the License, or
+## (at your option) any later version.
+##
+## This program is distributed in the hope that it will be useful,
+## but WITHOUT ANY WARRANTY; without even the implied warranty of
+## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+## GNU General Public License for more details.
+##
+## You should have received a copy of the GNU General Public License
+## along with this program; if not, write to the Free Software
+## Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
+
+## -*- texinfo -*-
+## @deftypefn {Function File} {@var{S} = } qtdecomp (@var{I})
+## @deftypefnx {Function File} {@var{S} = } qtdecomp (@var{I},@var{threshold})
+## @deftypefnx {Function File} {@var{S} = } qtdecomp (@var{I},@var{threshold},@var{mindim})
+## @deftypefnx {Function File} {@var{S} = } qtdecomp (@var{I},@var{threshold},@var{[mindim maxdim]})
+## @deftypefnx {Function File} {@var{S} = } qtdecomp (@var{I},@var{fun})
+## @deftypefnx {Function File} {@var{S} = } qtdecomp (@var{I},@var{fun},@var{P1},@var{P2},...)
+## Performs quadtree decomposition
+##
+## qtdecomp decomposes a square image @var{I} into four equal-sized
+## blocks. Then it performs some kind of test on each block to decide if
+## it should decompose them further. This process is repeated
+## iteratively until there's no block left to be decomposed.
+##
+## Note that blocks are not decomposed if their dimensions are not even.
+##
+## The output is a sparse matrix whose non-zero elements determine the
+## position of the block (the element is at top-left position in the
+## block) and size of each block (the value of the element determines
+## length of a side of the square-shaped block).
+##
+## S = qtdecomp(I) decomposes an intensity image @var{I} as described
+## above. By default it doesn't split a block if all elements are equal.
+##
+## S = qtdecomp(I, threshold) decomposes an image as decribed, but only
+## splits a block if the maximum value in the block minus the minimum
+## value is greater than @var{threshold}, which is a value between 0 and
+## 1. If @var{I} is of class uint8, @var{threshold} is multiplied by 255
+## before use. Also, if@var{I} is of class uint16, @var{threshold} is 
+## multiplied by 65535.
+##
+## S = qtdecomp(I, threshold, mindim) decomposes an image using the
+## @var{threshold} as just described, but doesn't produce blocks smaller
+## than mindim.
+##
+## S = qtdecomp(I, threshold, [mindim maxdim]) decomposes an image as
+## described, but produces blocks that can't be bigger than maxdim. It
+## decomposes to maxdim even if it isn't needed if only @var{threshold}
+## was considered.
+##
+## S = qtdecomp(I, fun) decomposes an image @var{I} and uses function
+## @var{fun} to decide if a block should be splitted or not. @var{fun}
+## is called with a m-by-m-by-k  array of m-by-m blocks to be
+## considered, and should return a vector of size k, whose elements
+## represent each block in the stacked array. @var{fun} sets the
+## corresponding value to 1 if the block should be split, and 0
+## otherwise.
+##
+## S = qtdecomp(I, fun, ...) behaves as qtdecomp(I, fun) but passes
+## extra parameters to @var{fun}.
+##
+## @end deftypefn
+
+## Author:  Josep Mones i Teixidor <jmones@puntbarra.com>
+
+function S = qtdecomp(I, p1, varargin)
+  if (nargin<1)
+    usage("S=qtdecomp(I)");
+  endif
+  
+  if (!ismatrix(I) || size(I,1)!=size(I,2))
+    error("qtdecomp: I should be a square matrix.");
+  endif
+
+  ## current size (assumed to be square)
+  curr_size=size(I,1);
+
+  ## initial mindim to a sensible value
+  mindim=1;
+ 
+  ## sensible default maxdim value
+  maxdim=curr_size;
+
+  if (nargin<2)
+    ## Initialize decision method variable
+    ## We could have implemented threshold as a function and use an
+    ## uniform interface (function handle) to decide whether to split or
+    ## not blocks. We have decided not to do so because block
+    ## rearrangement that is needed as a parameter to functions is
+    ## expensive and because it needs NDArray, which is only available
+    ## in Octave CVS.
+    decision_method=0;
+  elseif (isreal(p1))
+    ## p1 is threshold
+    threshold=p1;
+    decision_method=1;
+
+    if(strcmp(typeinfo(I), 'uint8 matrix'))
+      threshold*=255;
+    elseif(strcmp(typeinfo(I), 'uint16 matrix'))
+      threshold*=65535;
+    endif
+
+    if (nargin>3)
+      usage("S=qtdecomp(I,threshold,mindim),        \
+	  S=qtdecomp(I,threshold,[mindim maxdim])");
+    elseif (nargin==3)
+      dims=varargin{1};
+      if (isvector(dims)&&length(dims)==2)
+	mindim=dims(1);
+	maxdim=dims(2);
+      elseif (isreal(dims))
+	mindim=dims;
+      else
+	error("qtdecomp: third parameter must be 'mindim' or '[mindim maxdim]'");
+      endif
+      ## we won't check if mindim or maxdim are powers of 2. It's too
+      ## restrictive and don't need it at all.
+    endif
+    
+  elseif (strcmp(typeinfo(p1),"function handle") || \
+	  (isstr(p1) && strncmp(p1, "inline_func", 11)))
+    ## function handles seem to return true to isscalar
+    fun=p1;
+    decision_method=2;
+  else
+    error("qtdecomp: second parameter must be a integer (threshold) or a function handle (fun).");
+  endif
+  
+  ## initialize results matrices
+  res=[];
+
+  ## bool to flag end
+  finished=false;
+
+  ## array of offsets to blocks to evaluate
+  offsets=[1,1];
+
+  if (maxdim<mindim)
+    error("qtdecomp: mindim must be smaller than maxdim.");
+  endif
+
+  ## See if we have to split a minimum regarless other considerations.
+  if (maxdim<curr_size)
+    initial_splits=ceil(log2(curr_size/maxdim));
+    if(initial_splits>0)
+      divs=2^initial_splits;
+      if (rem(curr_size,divs)!=0)
+	error("qtdecomp: Can't decompose I enough times to fulfill maxdim requirement.");
+      endif
+      ## update curr_size
+      curr_size/=divs;
+      if(curr_size<mindim)
+	error("qtdecomp: maxdim restriction collides with mindim restriction.");
+      endif
+      els=transpose([0:divs-1]*curr_size+1);
+      offsets=[kron(els,ones(divs,1)), kron(ones(divs,1),els)];
+    endif
+  endif
+
+  while(!finished && rows(offsets)>0)
+    ## check other ending conditions:
+    ## is size is odd?
+    ## is splitted size < than mindim?
+    if ((rem(curr_size,2)!=0)||((curr_size/2)<mindim))
+      ## can't continue, lets add current evaluation blocks to results
+      res=[res; offsets, ones(size(offsets,1),1)*curr_size];
+      finished = true;
+    else
+      if (decision_method<2)
+	db=logical(ones(rows(offsets),1));
+	for r=1:rows(offsets)
+	  o=offsets(r,:);
+	  fo=offsets(r,:)+curr_size-1;
+
+	  if(decision_method==0)
+	    ## is everything equal?
+	    if (all(I(o(1),o(2))==I(o(1):fo(1),o(2):fo(2))))
+	      db(r)=0;
+	    endif
+	  else
+	    ## check threshold
+	    t=I(o(1):fo(1),o(2):fo(2));
+	    t=t(:);
+	    if ((max(t)-min(t))<=threshold)
+	      db(r)=0;
+	    endif
+	  endif
+	endfor
+      elseif(decision_method==2)
+	## function handle decision method
+	## build blocks
+	b=zeros(curr_size,curr_size,rows(offsets));
+	rbc=offsets(:,1:2)+curr_size-1;
+	for r=1:rows(offsets)
+	  b(:,:,r)=I(offsets(r,1):rbc(r,1),offsets(r,2):rbc(r,2));
+	endfor
+
+	db=feval(fun, b, all_va_args);
+      else
+	error("qtdecomp: execution shouldn't reach here. Please report this as a bug.");
+      endif
+
+      ## Add blocks that won't divide to results
+      nd=offsets(find(!db),:);
+      res=[res; nd, ones(size(nd,1),1)*curr_size];
+      
+      ## Update curr_size for next iteration
+      curr_size/=2;
+      
+      ## Prepare offsets for next iteration
+      otemp=offsets(find(db),:);
+      hs=ones(rows(otemp),1)*curr_size;
+      zs=zeros(size(hs));
+      offsets=[otemp;otemp+[hs,zs];otemp+[zs,hs];otemp+[hs,hs]];
+    endif
+  endwhile
+
+  S=sparse(res(:,1),res(:,2),res(:,3),size(I,1),size(I,2));
+endfunction
+
+
+%!demo
+%! full(qtdecomp(eye(8)))
+%! %It finds 2 big blocks of 0 and it decomposes further where 0 and 1 are mixed.
+
+
+%!# Test if odd-sized limits split
+%!assert(full(qtdecomp(eye(5))), reshape([5,zeros(1,24)],5,5));
+%!assert(full(qtdecomp(eye(6))), repmat(reshape([3,zeros(1,8)],3,3),2,2));
+
+%!shared A, B1, B2, B4, O2, O4, Req, R1, R10, R10m2
+%! A=[ 1, 4, 2, 5,54,55,61,62;
+%!     3, 6, 3, 1,58,53,67,65;
+%!     3, 6, 3, 1,58,53,67,65;
+%!     3, 6, 3, 1,58,53,67,65;
+%!    23,42,42,42,99,99,99,99;    
+%!    27,42,42,42,99,99,99,99;    
+%!    23,22,26,25,99,99,99,99;    
+%!    22,22,24,22,99,99,99,99];
+%! B1=[1];
+%! B2=[2,0;0,0];
+%! B4=reshape([4,zeros(1,15)],4,4);
+%! O2=ones(2,2);
+%! O4=ones(4,4);
+%! Req=[O4,O4; 
+%!      [O2,B2;O2,O2],B4];
+%! R1=[O4,O4;
+%!    [O2,B2;B2,O2],B4];
+%! R10=[[B4,[B2,B2;B2,B2]];
+%!      [[O2,B2;B2,B2],B4]];
+%! R10m2=[[B4,[B2,B2;B2,B2]];
+%!        [[B2,B2;B2,B2],B4]];
+
+%!# Test 'equal' method
+%!test
+%! a=ones(2,2);
+%! b=[2,0;0,0];
+%! assert(full(qtdecomp(eye(4))), [a,b;b,a]);
+
+%!assert(full(qtdecomp(A)), Req);
+
+%!# Test 'threshold' method
+%!assert(full(qtdecomp(A,0)), Req);
+%!assert(full(qtdecomp(A,1)), R1);
+%!assert(full(qtdecomp(A,10)), R10);
+%!assert(full(qtdecomp(A,10,2)), R10m2);
+%!assert(full(qtdecomp(A,100,[2, 4])), [B4,B4;B4,B4]);
+
+%!# Test 'fun' method
+%!# Can't define a function inside tests because Octave
+%!# complains about nested functions. If you know how to
+%!# do it please add function handle tests.
+
+%!# inline version
+%!# We use [size(A),1](3) as a trick since size(A,3) fails if A is 2D.
+
+%!# no params
+%!test
+%! first_eq=inline("(A(1,1,:)!=(54*ones(1,1,[size(A),1](3))))(:)","A");
+%! assert(full(qtdecomp(A,first_eq)),[O4,B4;O4,O4]); 
+
+%!# 1 param
+%!test
+%! first_eq=inline("(A(1,1,:)!=(c*ones(1,1,[size(A),1](3))))(:)","A","c");
+%! assert(full(qtdecomp(A,first_eq,54)),[O4,B4;O4,O4]); 
+
+%!# 3 params
+%!test
+%! first_eq=inline("(A(1,1,:)!=((c1+c2+c3)*ones(1,1,[size(A),1](3))))(:)","A","c1","c2","c3");
+%! assert(full(qtdecomp(A,first_eq,4,40,10)),[O4,B4;O4,O4]); 
+
+
+
+%
+% $Log$
+% Revision 1.1  2004/08/09 01:48:54  jmones
+% Added qtdecomp: quadtree decomposition
+%
+%
+