buildskelmapfast.sci

小波分解源代码

function [skellist,skelptr,skellen] = BuildSkelMapFast(maxmap)
// BuildSkelMapFast -- Chain together Ridges of Wavelet Transform
//  Usage
//    [skellist,skelptr,skellen] = BuildSkelMapFast(maxmap)
//  Inputs
//    maxmap    matrix MMWT
//  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
//    CWT, WTMM, PlotSkelMap, BuildSkelMap, ExtractRidge
//
//  Copyright Aldo I Maalouf

	[n,nscale,nvoice,noctave] = CalcCWTPars(size(maxmap));

	nchain = 0;
	chains = zeros(size(maxmap));
	count  = 0;

	chainset = [];
	currpos  = [];
	
	for k=1:nscale,		
		if k/10-floor(k/10)==1,
			mtlb_fprintf('BuildSkelMap: at scale %i:',k),
			mtlb_fprintf('nchains %i, #chainset %i\n',nchain,length2(chainset)),
		end
	
		// match existing features with finer scale features
		mindist = n/2 * 2^(-k/nvoice);
		
		if ~isempty(chainset),
			for l=1:length2(chainset),
				ipos = currpos(l);
				j = find(maxmap(:,k))'; 
				circdist   = min([ abs(j-ipos) ; abs(j-ipos+n); abs(j-ipos-n) ]);
				[dist,pos] = min(circdist);
	
				if isempty(pos),
					kill = 1;
				elseif dist(1) > mindist,
					kill = 1;
				else
					ipos = j(pos(1)); 
					kill = 0;
					maxmap(ipos,k)   = 0;
				end  // in principle could also test amplitude.
				
				if kill,  // kill current chain
					chains(chainset(l),k) = 0;
					chainset(l) = 0;
					currpos(l)  = 0;
				else      // grow current chain.	
					chains(chainset(l),k) = ipos;
					count = count+1;
					currpos(l) = ipos;
				end
			end
		end

		// look for roots of new chains, add them to chain list

    	j = find(maxmap(:,k));
		if ~isempty(j),
			ochain = nchain+1;
			for i=1:length2(j),
				nchain = nchain+1;
				chains(nchain,k) = j(i);
				count = count+1;
			end
	
			chainset = [ chainset (ochain:nchain) ];
			currpos  = [ currpos  chains(ochain:nchain,k)' ];
		end

		// clear out dead wood

		if ~isempty(chainset),
			j = find(chainset);
			chainset = chainset(j);
			currpos  = currpos(j) ;
		end

	end
	
	// built 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,

		if ch/10-floor(ch/10)==1,
			mtlb_fprintf('BuildSkelMap: packing chain %i\n',ch),
		end
	
		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