skelmap.sci

小波分解源代码

function [skellist,skelptr,skellen] = SkelMap(maxmap)
// SkelMap -- Chain together Ridges of Wavelet Transform
//  Usage
//    [skellist,skelptr,skellen] = SkelMap(maxmap)
//  Inputs
//    maxmap    matrix MM_RWT
//  Outputs
//    skellist  storage for list of chains
//    skelptr   vector of length nchain --pointers
//              to head of chain
//    skellen   vector of length nchain -- length of skellists
//
//  Description
//    A chain is a list of maxima at essentially the same position
//    across a range of scales. 
//    It is identified from the maxmap data structure output by WTMM
//    by finding a root at coarse scales and identifying the closest
//    maxima at the next finest scale.
//    NO PROVISION IS MADE FOR 'terminating' A CHAIN before the
//    finest scale is reached.
//
//    nchain = len(skellen) chains are found.
//    A chain data structure is a list of scale-location pairs
//    All chains are stored together in skellist.
//    The k-th list begins in skellist at skelptr(k)
//    The k-th list has length skellen(k)
//
//  See Also
//    RWT, MM_RWT, PlotSkelMap, ExtractRidge
//
//  Copyright Aldo I Maalouf

	[n,nscale] = size(maxmap);
	noctave = floor(log2(n))-5;
	nvoice  = nscale/noctave;

	nchain = 0;
	m=size(maxmap);
	chains = zeros(m(1),m(2));
	count  = 0;
		
	while any(any(maxmap)), //start new chain

	
	    [i,j] = find(maxmap);
	  
		iscale = j(1);
		ipos   = i(1);
		
		nchain = nchain+1;
		chains(nchain,iscale) = ipos;
		maxmap(ipos,iscale) = 0;
		count = count+1;
			
		while(iscale < nscale)  // pursue rest of chain
			iscale = iscale+1;
			j = find(maxmap(:,iscale))'; 
			circdist   = mtlb_min([ abs(j-ipos) ; abs(j-ipos+n); abs(j-ipos-n) ]);
			[dist,pos] = mtlb_min(circdist);
			if ~isempty(pos),
				ipos = j(pos(1)); 
				chains(nchain,iscale) = ipos;    
				
				maxmap(ipos,iscale)   = 0;
				count = count+1;
			else
				iscale = nscale;
			end
		end
	
	end
	
	// packed lists of chain structures

	rptr = zeros(1,n); 
	rlen = zeros(1,n);
	pchain = 1; qchain = 0;
	store = zeros(1,2*count);

	for ch =1:nchain,

		
	
		rptr(ch) = pchain;
		j = find(chains(ch,:));
		iscale = j(1);
		rlen(ch) = length2(j);
		ix  = iscale:(iscale+rlen(ch)-1);
		vec = [ ix ; chains(ch,ix)];
		qchain = pchain + (2*rlen(ch)-1);
		store(pchain:qchain) = vec(:)';  
	    pchain = qchain+1;
	end
	
	skelptr  = waverow(rptr(1:nchain));
	skellen  = waverow(rlen(1:nchain));
	skellist = waverow(store(1:qchain));
	endfunction