peaktransformdpimage.m

基于peak transform的高效图像编码。

if(0<1)
	fp=fopen('barb.512', 'rb');
	B=fscanf(fp, '%c', [512, 512]);
	fclose(fp);
	A=double(B');
end
[RTN, CTN]=size(A);

h0=[0.037828, -0.023849, -0.110624, 0.377402, 0.852699, 0.377402, -0.110624, -0.023849, 0.037828]  ;
h1=[0.064539, -0.040689,-0.418092, 0.788486, -0.418092, -0.040689, 0.064539];

	[ANA]=hvanalysis(h0, h1, A , 1,1, RTN, CTN, 1);
	[ANA]=hvanalysis(h0, h1, ANA, 1,1, RTN, CTN, -1);
	[ANA]=hvanalysis(h0, h1, ANA, 1,1, RTN/2, CTN/2, 1);
	[ANA]=hvanalysis(h0, h1, ANA, 1,1, RTN/2, CTN/2, -1);
	A = ANA(1:(RTN/4), 1:(CTN/4));


if(0 > 1)

	
	[ANA]=hvanalysis(h0, h1, ANA, 1,1, RTN/4, CTN/4, 1);




	[ANA]=hvanalysis(h0, h1, ANA, 1,1, RTN/4, CTN/4, -1);
end




[RTN, CTN]=size(A);
disp([RTN, CTN]);
pause;

%%----------------------parameter ------------------------------------------


LN = 300;
N = CTN;

QP = 8;

avgEnergyGain = QP/2;
peakThr = QP;


%%----------------------memory ------------------------------------------
orgF = zeros(1, CTN);

orgPeakMap = zeros(1, N/2);
winSize = 5;
numStates = 2^winSize;
upLinks = zeros(numStates, N) + 1;
ptrUpLinks = zeros(numStates, 1) + 1;
downLinks = zeros(numStates, N) + 1;
ptrDownLinks = zeros(numStates, 1) + 1;
baseValUpLinks = zeros(numStates, 1);
baseValDownLinks = zeros(numStates, 1);

G = orgF * 0;
orgFX = zeros(3, N);

flagUpDown = zeros(numStates, 1);
peakFlag = zeros(numStates, N/2);

posMap = 0:(numStates-1);
newPosMap = 0:(numStates-1);

testUpLink0 = orgF * 0;
testDownLink0 = orgF * 0;
testUpLink1 = orgF * 0;
testDownLink1 = orgF * 0;

oldUpLink0 = orgF * 0;
oldDownLink0 = orgF * 0;
oldUpLink1 = orgF * 0;
oldDownLink1 = orgF * 0;

peakMap = zeros(1, N/2);
truePeaks = zeros(1, N/2);

linkEnergy = zeros(1, numStates);


outputA = A * 0;
peakMapA = zeros(RTN, CTN/2);
PTStat = zeros(RTN, 6);

energyRatio = 1;
ct=1;

for LN=1:1:RTN, 

	%%---reset ---

		orgPeakMap = orgPeakMap * 0;

		upLinks = upLinks * 0;
		ptrUpLinks = ptrUpLinks * 0 + 1;
		downLinks = downLinks * 0;
		ptrDownLinks = ptrDownLinks * 0  + 1;

		G = orgF * 0;

		flagUpDown = flagUpDown * 0;
		peakFlag = peakFlag * 0;

		posMap = 0:(numStates-1);
		newPosMap = 0:(numStates-1);

		testUpLink0 = orgF * 0;
		testDownLink0 = orgF * 0;
		testUpLink1 = orgF * 0;
		testDownLink1 = orgF * 0;

		oldUpLink0 = orgF * 0;
		oldDownLink0 = orgF * 0;
		oldUpLink1 = orgF * 0;
		oldDownLink1 = orgF * 0;

		peakMap = peakMap * 0;
		truePeaks = truePeaks * 0;



	orgF = A(LN, :);

	%%----------------------select peaks ----------------------------------------
	%
	maxPeaks = N/2 * 0.6;

	[FANA]=hvanalysis(h0, h1, orgF, 1, 1, 1, N/1, 1);	

	numPeaks = 0;
	for k=(N/2+3):(N-3),
		if(abs(FANA(k)) > peakThr)
			pos = (k - N/2);
			orgPeakMap(pos) = 1;
		end
	end

	orgPeaks = zeros(1, N/2);
	numOrgPeaks = 0;

	for k=3:(N/2),
		if(orgPeakMap(k) > 0)
			numOrgPeaks = numOrgPeaks + 1;
			orgPeaks(numOrgPeaks) = k;
		end
	end

	numOrgPeaks = numOrgPeaks + 1;
	orgPeaks(numOrgPeaks) = N/2;

	if(numOrgPeaks > 10)

		%%----------------------start up -----------------------------
		%

		pos0 = 1;
		pos1 = 2 * orgPeaks(1);
		for ptrState=1:numStates,
			upLinks(ptrState, pos0:pos1) = orgF(pos0:pos1);
			ptrUpLinks(ptrState) = pos1;
			flagUpDown(ptrState) = 1;
			
		end


		for ptrStateVal=0:(numStates-1),

			ptrState = ptrStateVal + 1;

			for ptrOrgPeaks=1:winSize,

				pos0 = 2 * orgPeaks(ptrOrgPeaks);
				pos1 = 2 * orgPeaks(ptrOrgPeaks + 1);

				%%-- set the up / down flag --
				if(bitand(ptrStateVal, 2^(ptrOrgPeaks-1)) > 0)
					flagUpDown(ptrState) = (-1) * flagUpDown(ptrState);
					peakFlag(ptrState, ptrOrgPeaks) = 1;
				end


				%%-- copy the data --
				if(flagUpDown(ptrState) > 0)		
					currPos = ptrUpLinks(ptrState);
					upLinks(ptrState, currPos:(currPos + pos1 - pos0)) = orgF(pos0:pos1) - orgF(pos0) + upLinks(ptrState, currPos);
					ptrUpLinks(ptrState) = currPos + pos1 - pos0;
				else
					currPos = ptrDownLinks(ptrState);
					downLinks(ptrState, currPos:(currPos + pos1 - pos0)) = orgF(pos0:pos1) - orgF(pos0) + downLinks(ptrState, currPos);
					ptrDownLinks(ptrState) = currPos + pos1 - pos0;
				end

			end
		end


		%%----------------------Dynamic programming-----------------------------
		%


		for ptrOrgPeaks=(winSize+1):(numOrgPeaks-1),

			newPosMap = posMap * 0;

			for ptrStateVal=0:2:(numStates-1),
				ptrState = ptrStateVal + 1;

				linkIndex0 = posMap(ptrState) + 1;
				linkIndex1 = posMap(ptrState + 1) + 1;

				pos0 = 2 * orgPeaks(ptrOrgPeaks);
				pos1 = 2 * orgPeaks(ptrOrgPeaks + 1);

				flag0 = peakFlag(linkIndex0, :);
				flag1 = peakFlag(linkIndex1, :);

				oldUpLink0 = upLinks(linkIndex0, :);
				ptrOldUpLink0 = ptrUpLinks(linkIndex0);
				oldDownLink0 = downLinks(linkIndex0, :);
				ptrOldDownLink0 = ptrDownLinks(linkIndex0);

				oldUpLink1 = upLinks(linkIndex1, :);
				ptrOldUpLink1 = ptrUpLinks(linkIndex1);
				oldDownLink1 = downLinks(linkIndex1, :);
				ptrOldDownLink1 = ptrDownLinks(linkIndex1);

				oldFlagUpDown0 = flagUpDown(linkIndex0);
				oldFlagUpDown1 = flagUpDown(linkIndex1);


				for  k=0:1,
				
					if(k==0) 
						linkFlip = 1;
					else
						linkFlip = -1;
					end
					 
				%% the first link ----

					currFlag = oldFlagUpDown0 * linkFlip;

					if(currFlag > 0)
						currPos = ptrOldUpLink0;
						testUpLink0(:) = oldUpLink0;
						testUpLink0(currPos:(currPos + pos1 - pos0)) = orgF(pos0:pos1) - orgF(pos0) + testUpLink0(currPos);
						ptrTestUpLink0 = currPos + pos1 - pos0;

						testDownLink0(:) = oldDownLink0;
						ptrTestDownLink0 = ptrOldDownLink0;
					else
						testUpLink0(:) = oldUpLink0;
						ptrTestUpLink0 = ptrOldUpLink0;

						currPos = ptrOldDownLink0;
						testDownLink0(:) = oldDownLink0;
						testDownLink0(currPos:(currPos + pos1 - pos0)) = orgF(pos0:pos1) - orgF(pos0) + testDownLink0(currPos);
						ptrTestDownLink0 = currPos + pos1 - pos0;

					end


				%% the second link---

					currFlag = oldFlagUpDown1 * linkFlip;

					if(currFlag > 0)
						currPos = ptrOldUpLink1;
						testUpLink1(:) = oldUpLink1;
						testUpLink1(currPos:(currPos + pos1 - pos0)) = orgF(pos0:pos1) - orgF(pos0) + testUpLink1(currPos);
						ptrTestUpLink1 = currPos + pos1 - pos0;

						testDownLink1(:) = oldDownLink1;
						ptrTestDownLink1 = ptrOldDownLink1;
					else
						testUpLink1(:) = oldUpLink1;
						ptrTestUpLink1 = ptrOldUpLink1;

						currPos = ptrOldDownLink1;
						testDownLink1(:) = oldDownLink1;
						testDownLink1(currPos:(currPos + pos1 - pos0)) = orgF(pos0:pos1) - orgF(pos0) + testDownLink1(currPos);
						ptrTestDownLink1 = currPos + pos1 - pos0;

					end

			
				%%-- Find the energy ----
					
					%[En0] = findLinkEnergy(testUpLink0, ptrTestUpLink0, testDownLink0, ptrTestDownLink0);
					%[En1] = findLinkEnergy(testUpLink1, ptrTestUpLink1, testDownLink1, ptrTestDownLink1);

					En0 = 0;
					len = ptrTestUpLink0 + ptrTestDownLink0 - 1;
					sft = testUpLink0(ptrTestUpLink0) - testDownLink0(1);
					G(1:len) = [testUpLink0(1:ptrTestUpLink0), testDownLink0(2:ptrTestDownLink0) + sft];
					pa = 4;
					pb = len-3;
					if(pb >=pa)
						for mm = pa : 2 : pb,
							ac = G((mm-3):(mm+3)) * (h1');
							En0 = En0 + abs(ac);
						end
					end

					if(len ~= pos1)
						disp([ptrStateVal, k, 0]);
						pause;
					end

					En1 = 0;
					len = ptrTestUpLink1 + ptrTestDownLink1 - 1;
					sft = testUpLink1(ptrTestUpLink1) - testDownLink1(1);
					G(1:len) = [testUpLink1(1:ptrTestUpLink1), testDownLink1(2:ptrTestDownLink1) + sft];
					pa = 4;
					pb = len-3;
					if(pb >=pa)
						for mm = pa : 2 : pb,
							ac = G((mm-3):(mm+3)) * (h1');
							En1 = En1 + abs(ac);
						end
					end

					if(len ~= pos1)
						disp([ptrStateVal, k, 1]);
						pause;
					end

				%%-- choose the link----

					if (k==0)

						%	if(En0 < En1)
						if( (En0 - En1) < avgEnergyGain),

							%%-- choose link0 ---
							upLinks(linkIndex0, :) = testUpLink0;
							ptrUpLinks(linkIndex0) = ptrTestUpLink0; 
							
							downLinks(linkIndex0, :) = testDownLink0;
							ptrDownLinks(linkIndex0) = ptrTestDownLink0;

							flagUpDown(linkIndex0) = oldFlagUpDown0;
							linkEnergy(linkIndex0) = En0;

							peakFlag(linkIndex0, :) = flag0;
							peakFlag(linkIndex0, ptrOrgPeaks) = 0;

						else

							%%-- choose link1 ---
							upLinks(linkIndex0, :) = testUpLink1;
							ptrUpLinks(linkIndex0) = ptrTestUpLink1; 
							
							downLinks(linkIndex0, :) = testDownLink1;
							ptrDownLinks(linkIndex0) = ptrTestDownLink1;
							
							flagUpDown(linkIndex0) = oldFlagUpDown1;
							linkEnergy(linkIndex0) = En1;

							peakFlag(linkIndex0, :) = flag1;
							peakFlag(linkIndex0, ptrOrgPeaks) = 0;

						end

					%%--set the mapping---
						newIndex = floor(ptrStateVal / 2);
						newPosMap(newIndex + 1) = linkIndex0 - 1;
						

					else

						%	if(En0 < En1)
						if( (En0 - En1) < avgEnergyGain),

							%%-- choose link0 ---
							upLinks(linkIndex1, :) = testUpLink0;
							ptrUpLinks(linkIndex1) = ptrTestUpLink0; 

							downLinks(linkIndex1, :) = testDownLink0;
							ptrDownLinks(linkIndex1) = ptrTestDownLink0;

							flagUpDown(linkIndex1) = oldFlagUpDown0 * (-1);
							linkEnergy(linkIndex1) = En0;

							peakFlag(linkIndex1, :) = flag0;
							peakFlag(linkIndex1, ptrOrgPeaks) = 1;

						else

							%%-- choose link1 ---
							upLinks(linkIndex1, :) = testUpLink1;
							ptrUpLinks(linkIndex1) = ptrTestUpLink1; 
							
							downLinks(linkIndex1, :) = testDownLink1;
							ptrDownLinks(linkIndex1) = ptrTestDownLink1;

							flagUpDown(linkIndex1) = oldFlagUpDown1 * (-1);
							linkEnergy(linkIndex1) = En1;

							peakFlag(linkIndex1, :) = flag1;
							peakFlag(linkIndex1, ptrOrgPeaks) = 1;

						end

						%%--set the mapping---
						newIndex = floor(ptrStateVal / 2) + (numStates/2);
						newPosMap(newIndex + 1) = linkIndex1 - 1;

					end

				end		%%--end 0/1 ---

			end		%%-- for all links ---


			%%----------

			posMap = newPosMap;


			disp([LN, ptrOrgPeaks, numOrgPeaks, energyRatio, ct]);


		end


		%%----Finalize ------
			[minEnergy, minIndex] = min(linkEnergy);
			
			ptrTestUpLink = ptrUpLinks(minIndex);
			ptrTestDownLink = ptrDownLinks(minIndex);
			testUpLink = upLinks(minIndex, :);
			testDownLink = downLinks(minIndex, :);

			len = ptrTestUpLink + ptrTestDownLink - 1;
			sft = testUpLink(ptrTestUpLink) - testDownLink(1);
			G(1:len) = [testUpLink(1:ptrTestUpLink), testDownLink(2:ptrTestDownLink) + sft];
			G(len:N) = orgF(len:N) - orgF(len) + G(len);

			tmpPeakFlag = peakFlag(minIndex, 1:numOrgPeaks);

			ptrTruePeaks = 0;
			for ptrOrgPeaks = 1: numOrgPeaks,
				if(tmpPeakFlag(ptrOrgPeaks) > 0)
					pos = orgPeaks(ptrOrgPeaks);
					ptrTruePeaks = ptrTruePeaks + 1;
					truePeaks(ptrTruePeaks) = pos;
					peakMap(pos) = 1;
				end
			end

			ptrTruePeaks = ptrTruePeaks + 1;
			truePeaks(ptrTruePeaks) = N/2;

		%%--use these peaks to transform data ---
			numTruePeaks = ptrTruePeaks;

			orgFX(1, :) = orgF;
			orgFX(2, :) = 1:N; 

			H = orgFX * 0;
			upLinkH = H;
			downLinkH = H;
			
			ptrUpLinkH = 1;
			ptrDownLinkH = 1;

			pos1 = 2 * truePeaks(1);
			upLinkH(1, 1:pos1) = orgF(1, 1:pos1);
			upLinkH(2, 1:pos1) = 1:pos1;
			upLinkH(3, :) = 0;

			ptrUpLinkH = pos1;
			linkFlag = 1;

			for k=1:(numTruePeaks-1),
				linkFlag = (-1) * linkFlag;

				if(linkFlag > 0)
					currPos = ptrUpLinkH;
					pos0 = 2 * truePeaks(k);
					pos1 = 2 * truePeaks(k + 1);

					upLinkH(1, currPos:(currPos + pos1 - pos0)) = orgFX(1, pos0:pos1) - orgFX(1, pos0) + upLinkH(1, currPos);
					upLinkH(2, (currPos+1):(currPos + pos1 - pos0)) = (pos0+1) : pos1;
					upLinkH(3, (currPos+1):(currPos + pos1 - pos0)) =  - orgFX(1, pos0) + upLinkH(1, currPos);

					ptrUpLinkH = currPos + pos1 - pos0; 
				else
					currPos = ptrDownLinkH;
					pos0 = 2 * truePeaks(k);
					pos1 = 2 * truePeaks(k + 1);

					downLinkH(1, currPos:(currPos + pos1 - pos0)) = orgFX(1, pos0:pos1) - orgFX(1, pos0) + downLinkH(1, currPos);
					downLinkH(2, (currPos+1):(currPos + pos1 - pos0)) = (pos0+1) : pos1;
					downLinkH(3, (currPos+1):(currPos + pos1 - pos0)) =  - orgFX(1, pos0) + downLinkH(1, currPos);

					ptrDownLinkH = currPos + pos1 - pos0; 
				end

			end

			len = ptrUpLinkH + ptrDownLinkH - 1 ;
			sft = upLinkH(1, ptrUpLinkH)  - downLinkH(1, 1);
			H(1, 1:len) = [upLinkH(1, 1:ptrUpLinkH), downLinkH(1, 2:ptrDownLinkH) + sft];
			H(2, 1:len) = [upLinkH(2, 1:ptrUpLinkH), downLinkH(2, 2:ptrDownLinkH)];
			H(3, 1:len) = [upLinkH(3, 1:ptrUpLinkH), downLinkH(3, 2:ptrDownLinkH) + sft];

			[GANA]=hvanalysis(h0, h1, G, 1, 1, 1, N/1, 1);	
			[HANA]=hvanalysis(h0, h1, H(1, :), 1, 1, 1, N/1, 1);	

			EG = sum(abs(GANA((N/2+1):N)));
			EH = sum(abs(HANA((N/2+1):N)));
			EF = sum(abs(FANA((N/2+1):N)));


			%%--inverse PT---
			PTFANA = zeros(1, N);

			for k=1:(N/2),
				pos = H(2, 2*k-1);
				x = HANA(k) - H(3, 2*k-1) * 1.412;
				ind = ceil(pos / 2);
				PTFANA(ind) = x;

				pos = H(2, 2*k);
				x = HANA(k+N/2);
				ind = ceil(pos / 2);
				PTFANA(ind+N/2) = x;
			end

			PTFANA = HANA;

		
		
			
			%%-- find the non-zeros;
			
			QF = round( FANA((N/2+1):N) / QP);
			NZF = abs(sign(QF));

			QH = round( PTFANA((N/2+1):N) / QP);
			NZH = abs(sign(QH));

			energyRatio = EH/EF
			ct = sum(peakMap);

			disp([sum(abs(QF)), sum(abs(QH)), sum(peakMap), floor(max(abs(G-H(1, :))))]);
			disp('----------------------------------------------');
			outputA(LN, :) = PTFANA;
			peakMapA(LN, :) = peakMap;
			PTStat(LN, 1) = EF;
			PTStat(LN, 2) = EG;
			PTStat(LN, 3) = EH;
			PTStat(LN, 4) = sum(abs(QF));
			PTStat(LN, 5) = sum(abs(QH));
			PTStat(LN, 6) = sum(peakMap);
		else
			disp(numOrgPeaks);
			outputA(LN, :) = FANA;
			peakMapA(LN, :) = 0;
			PTStat(LN, :) = 0;
		end
		
end

%%----------------------display ----------------------------------------
%

figure;
imagesc(outputA);
colormap(gray);


figure;
imagesc(peakMapA);
colormap(gray);



if(0>1)

	figure;
	plot(orgF, 'b.-');
	hold on;
	plot(G, 'r--');

	ym0 = max(abs(FANA((N/2+2):(N-2))));
	ym1 = max(abs(FANA((N/2+2):(N-2))));
	ym = max(ym0, ym1);

	figure;
	subplot(2, 1, 1);
	T = (2):(N/2-2);
	T = 2*T;
	plot(T, FANA((N/2+2):(N-2)), '.-');
	axis([1, 512, -ym, ym]);
	subplot(2, 1, 2);
	plot(T, GANA((N/2+2):(N-2)), '.-');
	axis([1, 512, -ym, ym]);


	figure;
	subplot(2, 1, 1);
	imagesc(peakFlag(:, 1:numOrgPeaks));
	colormap(gray);
	subplot(2, 1, 2);
	stem(peakMap, '.');

	disp(sum(peakMap));

	figure;
	plot(G, '.-');
	hold on;
	plot(H(1, :), 'rd-');

	figure;
	subplot(2, 1, 1);
	plot(FANA(1:N), '.-');
	subplot(2, 1, 2);
	plot(PTFANA(1:N), '.-');

end

AH2 = outputA;
MH2 = peakMapA;
SH2 = PTStat;

save AH2_data AH2;
save MH2_data MH2;
save SH2_data SH2;