%function expand (by interpolation), (somewhat) inverse of reduce

function gl1=expand(gl0,w)

[C,R]=size(gl0);
gl1=zeros(2*C-1,2*R-1);

%interpolation
gl1(2:2:2*C-2,2:2:2*R-2)=4*(gl0(1:C-1,1:R-1)*w(4,4)+gl0(2:C,1:R-1)*w(2,4)+gl0(1:C-1,2:R)*w(4,2)+gl0(2:C,2:R)*w(2,2));
gl1(3:2:2*C-3,2:2:2*R-2)=4*(gl0(1:C-2,1:R-1)*w(5,4)+gl0(2:C-1,1:R-1)*w(3,4)+gl0(3:C,1:R-1)*w(1,4)+gl0(1:C-2,2:R)*w(5,2)+gl0(2:C-1,2:R)*w(3,2)+gl0(3:C,2:R)*w(1,2));
gl1(2:2:2*C-2,3:2:2*R-3)=4*(gl0(1:C-1,1:R-2)*w(4,5)+gl0(1:C-1,2:R-1)*w(4,3)+gl0(1:C-1,3:R)*w(4,1)+gl0(2:C,1:R-2)*w(2,5)+gl0(2:C,2:R-1)*w(2,3)+gl0(2:C,3:R)*w(2,1));
gl1(3:2:2*C-3,3:2:2*R-3)=4*(gl0(1:C-2,1:R-2)*w(5,5)+gl0(2:C-1,1:R-2)*w(3,5)+gl0(3:C,1:R-2)*w(1,5)+gl0(1:C-2,2:R-1)*w(5,3)+gl0(2:C-1,2:R-1)*w(3,3)+gl0(3:C,2:R-1)*w(1,3)+gl0(1:C-2,3:R)*w(5,1)+gl0(2:C-1,3:R)*w(3,1)+gl0(3:C,3:R)*w(1,1));

%compute edges
temp=w(5,3)+w(5,1)+w(3,3)+w(3,1)+w(1,3)+w(1,1);
gl1(3:2:2*C-3,1)=(gl0(1:C-2,1)*w(5,3)+gl0(1:C-2,2)*w(5,1)+gl0(2:C-1,1)*w(3,3)+gl0(2:C-1,2)*w(3,1)+gl0(3:C,1)*w(1,3)+gl0(3:C,2)*w(1,1))/temp;
temp=w(5,3)+w(5,5)+w(3,3)+w(3,5)+w(1,3)+w(1,5);
gl1(3:2:2*C-3,2*R-1)=(gl0(1:C-2,R)*w(5,3)+gl0(1:C-2,R-1)*w(5,5)+gl0(2:C-1,R)*w(3,3)+gl0(2:C-1,R-1)*w(3,5)+gl0(3:C,R)*w(1,3)+gl0(3:C,R-1)*w(1,5))/temp;
temp=w(4,3)+w(4,1)+w(2,3)+w(2,1);
gl1(2:2:2*C-2,1)=(gl0(1:C-1,1)*w(4,3)+gl0(1:C-1,2)*w(4,1)+gl0(2:C,1)*w(2,3)+gl0(2:C,2)*w(2,1))/temp;
temp=w(4,3)+w(4,5)+w(2,3)+w(2,5);
gl1(2:2:2*C-2,2*R-1)=(gl0(1:C-1,R)*w(4,3)+gl0(1:C-1,R-1)*w(4,5)+gl0(2:C,R)*w(2,3)+gl0(2:C,R-1)*w(2,5))/temp;

temp=w(3,5)+w(1,5)+w(3,3)+w(1,3)+w(3,1)+w(1,1);
gl1(1,3:2:2*R-3)=(gl0(1,1:R-2)*w(3,5)+gl0(2,1:R-2)*w(1,5)+gl0(1,2:R-1)*w(3,3)+gl0(2,2:R-1)*w(1,3)+gl0(1,3:R)*w(3,1)+gl0(2,3:R)*w(1,1))/temp;
temp=w(3,5)+w(5,5)+w(3,3)+w(5,3)+w(3,1)+w(5,1);
gl1(2*C-1,3:2:2*R-3)=(gl0(C,1:R-2)*w(3,5)+gl0(C-1,1:R-2)*w(5,5)+gl0(C,2:R-1)*w(3,3)+gl0(C-1,2:R-1)*w(5,3)+gl0(C,3:R)*w(3,1)+gl0(C-1,3:R)*w(5,1))/temp;
temp=w(3,4)+w(1,4)+w(3,2)+w(1,2);
gl1(1,2:2:2*R-2)=(gl0(1,1:R-1)*w(3,4)+gl0(2,1:R-1)*w(1,4)+gl0(1,2:R)*w(3,2)+gl0(2,2:R)*w(1,2))/temp;
temp=w(3,4)+w(5,4)+w(3,2)+w(5,2);
gl1(2*C-1,2:2:2*R-2)=(gl0(C,1:R-1)*w(3,4)+gl0(C-1,1:R-1)*w(5,4)+gl0(C,2:R)*w(3,2)+gl0(C-1,2:R)*w(5,2))/temp;

%compute corners
temp=w(3,3)+w(3,1)+w(1,3)+w(1,1);
gl1(1,1)=(gl0(1,1)*w(3,3)+gl0(1,2)*w(3,1)+gl0(2,1)*w(1,3)+gl0(2,2)*w(1,1))/temp;
temp=w(3,3)+w(5,3)+w(3,1)+w(5,1);
gl1(2*C-1,1)=(gl0(C,1)*w(3,3)+gl0(C-1,1)*w(5,3)+gl0(C,2)*w(3,1)+gl0(C-1,2)*w(5,1))/temp;
temp=w(3,3)+w(3,5)+w(1,3)+w(1,5);
gl1(1,2*R-1)=(gl0(1,R)*w(3,3)+gl0(1,R-1)*w(3,5)+gl0(2,R)*w(1,3)+gl0(2,R-1)*w(1,5))/temp;
temp=w(3,3)+w(5,3)+w(3,5)+w(5,5);
gl1(2*C-1,2*R-1)=(gl0(C,R)*w(3,3)+gl0(C-1,R)*w(5,3)+gl0(C,R-1)*w(3,5)+gl0(C-1,R-1)*w(5,5))/temp;


