Example.m

clear variables;clc;

% open an image in Zeiss LSM format that is LZW encoded.
fid=fopen('LZW.lsm');
imageHeader=fread(fid,95950,'int8');
imageDataLZW=fread(fid,'int8'); fclose(fid); clear fid ans;

% decompress the LZW encoded stream (must be column vector of doubles)
tic;imageData=lzwDecompressM(imageDataLZW);tmp=toc;
display(['Decoding time for 512x512 image was ' num2str(tmp) ' seconds']); clear tmp;

% bring image data to uint8 format
imageData=double(reshape(imageData,[512 512])');
imageBuffer=zeros(512);

imageBuffer(:,1)=imageData(:,1);
for j=2:512
    imageBuffer(:,j)=imageBuffer(:,j-1)+imageData(:,j);
end
while min(min(imageBuffer))<0
    imageBuffer=imageBuffer+256.*(imageBuffer<0);
end
imageBuffer=uint32(imageBuffer);
imageBuffer=bitshift(bitshift(imageBuffer,24),-24);
imageData=uint8(imageBuffer); clear imageBuffer j;

% apply some filtering
imageFilteredData=imageData;
h1=logical(repmat([0 1 0],512,1));
h2=logical(repmat([0 0 1 0 0],512,1));
h3=logical(repmat([0 1 1 0],512,1));
h4=logical(repmat([0 0 1 1 0 0],512,1));


for j=1:512-size(h1,2)
    imageFilteredData(find(sum(squeeze(imageFilteredData(:,j:j+size(h1,2)-1)>0).*h1,2)==sum(h1,2) & sum(imageFilteredData(:,j:j+size(h1,2)-1)>0,2)==sum(h1,2)),j:j+size(h1,2)-1)=0;
end

for j=1:512-size(h2,2)
    imageFilteredData(find(sum(squeeze(imageFilteredData(:,j:j+size(h2,2)-1)>0).*h2,2)==sum(h2,2) & sum(imageFilteredData(:,j:j+size(h2,2)-1)>0,2)==sum(h2,2)),j:j+size(h2,2)-1)=0;
end

for j=1:512-size(h3,2)
    imageFilteredData(find(sum(squeeze(imageFilteredData(:,j:j+size(h3,2)-1)>0).*h3,2)==sum(h3,2) & sum(imageFilteredData(:,j:j+size(h3,2)-1)>0,2)==sum(h3,2)),j:j+size(h3,2)-1)=0;
end

for j=1:512-size(h4,2)
    imageFilteredData(find(sum(squeeze(imageFilteredData(:,j:j+size(h4,2)-1)>0).*h4,2)==sum(h4,2) & sum(imageFilteredData(:,j:j+size(h4,2)-1)>0,2)==sum(h4,2)),j:j+size(h4,2)-1)=0;
end


clear j h1 h2 h3 h4;

% bring the image data to compression formatting by horizontal
% differencing.
imageFilteredData=int16(imageFilteredData);
imageBuffer=zeros(512,'int16');

imageBuffer(:,1)=imageFilteredData(:,1);
imageBuffer(:,2:end)=imageFilteredData(:,2:end)-imageFilteredData(:,1:end-1);
imageFilteredData=uint8(imageFilteredData);

imageBuffer=imageBuffer+256.*int16(imageBuffer<0);
imageBuffer=uint16(imageBuffer);
imageBuffer=bitshift(bitshift(imageBuffer,8),-8);
imageBuffer=int16(imageBuffer);

imageBuffer=imageBuffer-256.*int16(imageBuffer>127);
imageBuffer=reshape(imageBuffer',[262144 1]);

tic;imageFilteredDataLZW=lzwCompressM(double(imageBuffer));tmp=toc;
display(['Encoding time for 512x512 image was ' num2str(tmp) ' seconds']); clear tmp imageBuffer;

fid=fopen('lzwFiltered.lsm','wb');
fwrite(fid,imageHeader,'int8');
fwrite(fid,imageFilteredDataLZW,'int8');
fclose(fid); clear fid ans imageHeader;