xx=(1:128);
[X Y]=meshgrid(xx,xx);
x=avsread('phtm.fld')';
x(find(x==1000 | x==5600))=0; 		%desk
x(find(x==26000))=0.096;		%heart
x(find(x==22000))=0.096;		%body
x(find(x==32000))=0.17;			%bone
x(find(x==12000))=0.022;		%lung
x=x+0.002*(sin(X/2)+sin(Y/2)).*(x==0.022 & X<64);
avswrite('phtmt.fld',x');

m=(Y>5 & Y<123 & X>20 & X<110);
avswrite('mask.fld',m);

% blurring filter sigma=0.7 or 0.8
