function Mats_02

m=4;
n=5;
M=MatConfig(m,n);

%filter for vertical reflections
Vref=MatSymV(M,m,n); 
V=(Vref==1:size(M,2));
V=V(Vref>=0);
M=M(:,Vref>=0);

%filter for horizontal reflections 
Href=MatSymH(M,m,n);
H=(Href==1:size(M,2));
H=H(Href>=0); 
V=V(Href>=0);
M=M(:,Href>=0);

%filter for 180 degree rotations
Rref=MatSymR(M,m,n);
R180=(Rref==1:size(M,2));
R180=R180(Rref>=0);
H=      H(Rref>=0); 
V=      V(Rref>=0);
M=      M(:,Rref>=0);

if m==n
%filter for 90 degree rotations
Rref=MatSymR2(M,m,n);
R90=(Rref==1:size(M,2));
R90=  R90(Rref>=0);
R180=R180(Rref>=0);
H=      H(Rref>=0); 
V=      V(Rref>=0);
M=      M(:,Rref>=0);
%Now remove configs that differ by a reflection after a 90 degree rotation.
Rref=logical(MatSymR2b(M,m,n))
R90=  R90(Rref);
R180=R180(Rref);
H=      H(Rref); 
V=      V(Rref);
M=      M(:,Rref);
[1:length(V); R90; R180; H; V]
else
[1:length(V); R180; H; V]
end

S=ceil(size(M,2)/36);
count=0;
for ks=1:S
    figure('color',[1 1 1])
    for k=1:36
        count=count+1;
        subplot(6,6,k)
        MatPlot(reshape(M(:,count),m,n));
        title(num2str(count))
        axis equal
        axis off
        if count==size(M,2)
            break
        end
    end
end

%--------------------------------------------------------------------------
function Vref=MatSymV(M,m,n)
%find configs symmetric by reflection about the vertical

%keep track of what's already been paired off
N=size(M,2);
Vref=zeros(1,N);    Iv=true(1,N);
for k=1:N
    %vertical reflections
    if Iv(k) %k hasn't been paired yet
        Ck=reshape(M(:,k),m,n);
        Ck=flipud(Ck);
        I=find(Ck==1);
        for k3=I
            Ck(k3-1)=1;
            Ck(k3)=-1;
        end
        for k2=k:N
            if Iv(k2)   %k2 hasn't been paired yet
                if all(Ck(:)==M(:,k2))
                    Vref(k)=-1; Iv(k)= false;
                    Vref(k2)=k; Iv(k2)=false;
                    break
                end
            end
        end
    end
end

%--------------------------------------------------------------------------
function Href=MatSymH(M,m,n)
%find configs symmetric by reflection about the horizontal

%keep track of what's already been paired off
N=size(M,2);
Href=zeros(1,N);    
Ih=true(1,N);
for k=1:N
    if Ih(k) %k hasn't been paired yet
        Ck=reshape(M(:,k),m,n);
        Ck=fliplr(Ck);
        I=find(Ck==2);
        for k3=I
            Ck(k3-m)=2;
            Ck(k3)=-1;
        end
        for k2=k:N
            if Ih(k2)   %k2 hasn't been paired yet
                if all(all(Ck(:)==M(:,k2)))
                    Href(k)=-1; Ih(k)= false;
                    Href(k2)=k; Ih(k2)=false;
                    break
                end
            end
        end
    end
end

%--------------------------------------------------------------------------
function R180=MatSymR(M,m,n)
%find configs symmetric by a 180 degree rotation

%keep track of what's already been paired off
N=size(M,2);
R180=zeros(1,N);    
Ir=true(1,N);
for k=1:N
    if Ir(k) %k hasn't been paired yet
        Ck=reshape(M(:,k),m,n);
        Ck=rot90(rot90(Ck));
        I=find(Ck==2);
        for k3=I
            Ck(k3-m)=2;
            Ck(k3)=-1;
        end
        I=find(Ck==1);
        for k3=I
            Ck(k3-1)=1;
            Ck(k3)=-1;
        end
        for k2=k:N
            if Ir(k2)   %k2 hasn't been paired yet
                if all(all(Ck(:)==M(:,k2)))
                    R180(k)=-1; Ir(k)= false;
                    R180(k2)=k; Ir(k2)=false;
                    break
                end
            end
        end
    end
end

%--------------------------------------------------------------------------
function R90=MatSymR2(M,m,n)
%find configs symmetric by a 90 degree rotation

%keep track of what's already been paired off
N=size(M,2);
R90=zeros(1,N);    
Ir=true(1,N);
for k=1:N
    if Ir(k) %k hasn't been paired yet
        Ck=reshape(M(:,k),m,n);
        Ck=rot90(Ck);
        I1=find(Ck==1);
        I2=find(Ck==2);
        for k3=I1
            Ck(k3)=2;
        end
        for k3=I2
            Ck(k3-1)=1;
            Ck(k3)=-1;
        end
        for k2=k:N
            if Ir(k2)   %k2 hasn't been paired yet
                if all(all(Ck(:)==M(:,k2)))
                    R90(k)=-1; Ir(k)= false;
                    R90(k2)=k; Ir(k2)=false;
                    break
                end
            end
        end
    end
end

%--------------------------------------------------------------------------
function R90=MatSymR2b(M,m,n)
%M(:,R90) is a set of configs no two of which are related by a 90 degree
%rotation followed by a reflection.

%keep track of what's already been paired off
N=size(M,2);
R90=ones(1,N);    
Ir=true(1,N);
for k=1:N
    if Ir(k) %k hasn't been eliminated yet
        Ck1=reshape(M(:,k),m,n);
        Ck1=flipud(rot90(Ck1));
        I1=find(Ck1==1);
        I2=find(Ck1==2);
        for k3=I1
            Ck1(k3)=2;
        end
        for k3=I2
            Ck1(k3)=1;
        end
        %
        Ck2=reshape(M(:,k),m,n);
        Ck2=fliplr(rot90(Ck2));
        I1=find(Ck2==1);
        I2=find(Ck2==2);
        for k3=I1
            Ck2(k3)=-1;
            Ck2(k3-m)=2;
        end
        for k3=I2
            Ck2(k3)=-1;
            Ck2(k3-1)=1;
        end
        %
        for k2=k+1:N
            if Ir(k2)   %k2 hasn't been paired yet
                if all(all(Ck1(:)==M(:,k2))) || all(all(Ck2(:)==M(:,k2))) 
                    R90(k2)=0; Ir(k2)=false;
                    break
                end
            end
        end
    end
end

%--------------------------------------------------------------------------
function M=MatConfig(m,n)
%return raw configs.  no identifications because of reflection or rotation.
%
%   m  number of rows
%   n  number of columns

if floor(0.5*m*n)~=(0.5*m*n) %either m or n must be even.
    M=[];
    return
end
M=zeros(m*n,2^(0.5*m*n));
count=1;
for k=1:2^(0.5*m*n)
    C=zeros(m,n);
    flag=0;
    for kb=1:0.5*m*n
        I=find(C==0,1);
        c=ceil(I/m);
        r=I-m*(c-1);
        if bitget(k,kb) %go down
            if r==m || C(r+1,c)~=0, flag=1; break, end
            C(r,c)=1;
            C(r+1,c)=-1;
        else            %go right
            if c==n, flag=1; break, end
            C(r,c)=2;
            C(r,c+1)=-1;
        end
    end
    if flag, continue, end
    M(:,count)=C(:);
    count=count+1;
end
M=M(:,1:count-1);

%--------------------------------------------------------------------------
function MatPlot(C)
%create a plot of the configuration of mats corresponding to C.
[m,n]=size(C);
for k=1:m*n
    if C(k)==-1; continue, end
    c=ceil(k/m);
    r=k-m*(c-1);
    if C(k)==1
        line('xdata',c+[0 0 0 1],'ydata',[-2 -1 0 0]+m-r)
    else
        line('xdata',c+[0 0 1 2],'ydata',[-1 0 0 0]+m-r)
    end
end
line('xdata',[1 n+1 n+1],'ydata',[0 0 m]-1,'color',[0 0 0])
axis off
axis equal