%function PatternSequenceAlbedo()
% Hassebrook copyright 2007
%read in pattern sequence
clear all;
Nf = 2; % number of different frequencies
Nk(1)=1; % kc0=1 (k=0) first of Nf frequences
Np(1)=32; % number pattern shifts for first frequency
Nk(2)=16; % kc1=16 (k=1) second of Nf frequences
Np(2)=32; % number pattern shifts for second frequency
Pathname='Data3D' % path or folder name in reference to your default path
Filename='pattern' % name of files not including the index or suffix
Filesuffix='bmp' % suffix or image type
index=0; % first index is 0
Fullname=sprintf('%s%c%s%d%c%s',Pathname,'\',Filename,index,'.',Filesuffix)
A_bmp=double(imread(Fullname)); % load pattern#.bmp
[My, Nx, Pz] =size(A_bmp);
My
Nx
Pz
Albedo=zeros(My,Nx);
Phase0=zeros(My,Nx);
Phase1=zeros(My,Nx);
for k=1:Nf;
    Iq=zeros(My,Nx);
    Ii=zeros(My,Nx);
    for p=1:Np(k);
        Fullname=sprintf('%s%c%s%d%c%s',Pathname,'\',Filename,index,'.',Filesuffix)
        A_bmp=double(imread(Fullname)); % load pattern#.bmp
        Ar=A_bmp(:,:,1);
        % average in albedo
        % determin the sin and cos phase shift coefficients
        sincoef=sin(2*pi*(p-1)/Np(k));
        coscoef=cos(2*pi*(p-1)/Np(k));
        % average in weighted quadrature and inphase images
        % increment to next pattern
        index=index+1;
    end; % loop over p
    % divide quadrature and phase component by their number of shifts 
    if Np(k)>1
    end;
    % process first set of patterns into phase image
    if k==1
    end;
    % process second set of patterns into phase image 
    if k==2
    end
end; % loop over k
index
Albedo=Albedo/index;
figure(1);
imagesc(Albedo);
colormap gray;
figure(2);
imagesc(Phase0);
colormap gray;
figure(3);
imagesc(Phase1);
colormap gray;
P0(1:My)=Phase0(1:My,Nx/2);
P1(1:My)=Phase1(1:My,Nx/2);
m=1:My;
figure(4);
plot(m,P0,m,P1);
print -dtiff 'PhaseCrossSection.tif';
figure(5);
plot(m,P0/Nk(1),m,P1/Nk(2));
print -dtiff 'PhaseNormCrossSection.tif';
%
Albedo8=uint8(Albedo);
ImageI=Albedo8;
max(max(Albedo8))
min(min(Albedo8))
B_bmp(:,:,1)=Albedo8;
B_bmp(:,:,2)=Albedo8;
B_bmp(:,:,3)=Albedo8;
bmpimageC=B_bmp;
bmpimageI=B_bmp;
imwrite(B_bmp,'Albedo.bmp','BMP');
%
Albedo8=uint8(Phase0*255/(2*pi));
max(max(Albedo8))
min(min(Albedo8))
B_bmp(:,:,1)=Albedo8;
B_bmp(:,:,2)=Albedo8;
B_bmp(:,:,3)=Albedo8;
imwrite(B_bmp,'Phase0.bmp','BMP');
%
Albedo8=uint8(Phase1*255/(2*pi));
max(max(Albedo8))
min(min(Albedo8))
B_bmp(:,:,1)=Albedo8;
B_bmp(:,:,2)=Albedo8;
B_bmp(:,:,3)=Albedo8;
imwrite(B_bmp,'Phase1.bmp','BMP');
%
Albedo8=uint8(Phase00*255/(2*pi));
max(max(Albedo8))
min(min(Albedo8))
B_bmp(:,:,1)=Albedo8;
B_bmp(:,:,2)=Albedo8;
B_bmp(:,:,3)=Albedo8;
imwrite(B_bmp,'PhaseYP.bmp','BMP');

result=CPleastsqmethod();
%
out=getWorldCoordsB(ImageI,Phase00,32);
xw=out(:,:,1);
yw=out(:,:,2);
zw=out(:,:,3);
mat5file='fred';
result = mat5write(mat5file,bmpimageC,bmpimageI,xw,yw,zw);