function g = grayWaveletCompress(grayImage,WaveletOrder,scale,ratio)
% works for 512x512 images but can be adapted
[m,n] = size(grayImage);
% grayscale wavelet Mandrill compression
load(['Daub' num2str(WaveletOrder)]);
grayImage=double(grayImage);
grayImage.s = grayImage;grayImage.rd=0;grayImage.cd=0;
h.s=h;h.d=0;
g.s=g;g.d=0;
rh.s=rh;rh.d=WaveletOrder-1;
rg.s=rg;rg.d=WaveletOrder-1;
numEl=round(512*512/ratio);
WT = WaveTransform2D(grayImage,h,g,scale);
t = thresholdValue2D(WT,numEl);
WT = threshold2D(WT,t);
v = InvWaveTransform2D(WT,rh,rg,scale);
% trim columns at left that comes from the rd field and adjust length
g = v.s(:,v.rd+1:v.rd+n);
% trim rows
g = g(v.cd+1:v.cd+m,:);
g = uint8(g);
% works for 512x512 images but can be adapted
[m,n] = size(grayImage);
% grayscale wavelet Mandrill compression
load(['Daub' num2str(WaveletOrder)]);
grayImage=double(grayImage);
grayImage.s = grayImage;grayImage.rd=0;grayImage.cd=0;
h.s=h;h.d=0;
g.s=g;g.d=0;
rh.s=rh;rh.d=WaveletOrder-1;
rg.s=rg;rg.d=WaveletOrder-1;
numEl=round(512*512/ratio);
WT = WaveTransform2D(grayImage,h,g,scale);
t = thresholdValue2D(WT,numEl);
WT = threshold2D(WT,t);
v = InvWaveTransform2D(WT,rh,rg,scale);
% trim columns at left that comes from the rd field and adjust length
g = v.s(:,v.rd+1:v.rd+n);
% trim rows
g = g(v.cd+1:v.cd+m,:);
g = uint8(g);