clear
clc

nPhotos=3;
imgl1=imread('DSCN1037.jpg');
imgl1=imresize(imgl1,0.01);
%imgl1=im2double(imgl1);
imgl2=imread('DSCN1036.jpg');
imgl2=imresize(imgl2,0.01);
%imgl2=im2double(imgl2);
imgl3=imread('DSCN1035.jpg');
imgl3=imresize(imgl3,0.01);
%imgl3=im2double(imgl3);

[RowImgl ColImgl RGBImgl]=size(imgl1);

nPixels=RowImgl*ColImgl;

gZ1_R=reshape(imgl1(:,:,1),RowImgl*ColImgl,1);
gZ2_R=reshape(imgl2(:,:,1),RowImgl*ColImgl,1);
gZ3_R=reshape(imgl3(:,:,1),RowImgl*ColImgl,1);
gZ_R=[gZ1_R;gZ2_R;gZ3_R];

%coeffMtrx=sparse(nPixels*nPhotos+254,nPixels+256);
lnT_gz=[ones(nPixels,1).*log(0.5);ones(nPixels,1).*log(0.2);ones(nPixels,1).*log(0.02);zeros(254,1)];
RowIdx=(zeros(nPhotos*nPixels*2+254*3,1));
ColIdx=(zeros(nPhotos*nPixels*2+254*3,1));
val=zeros(nPhotos*nPixels*2+254*3,1);

for i=1:nPhotos
    for j=1:nPixels
        RowIdx(((i-1)*nPixels+j)*2-1,1)=(i-1)*nPixels+j;
        ColIdx(((i-1)*nPixels+j)*2-1,1)=gZ_R((i-1)*nPixels+j,1)+1;
        val(((i-1)*nPixels+j)*2-1,1)=1;
        
        RowIdx(((i-1)*nPixels+j)*2,1)=(i-1)*nPixels+j;
        ColIdx(((i-1)*nPixels+j)*2,1)=256+j;
        val(((i-1)*nPixels+j)*2,1)=-1;
      %  coeffMtrx((i-1)*nPixels+j,gZ_R((i-1)*nPixels+j,1))=1;
      %  coeffMtrx((i-1)*nPixels+j,256+j)=-1;
    end
end

for i=1:254
    RowIdx((nPhotos*nPixels)*2+i*3-2,1)=nPixels*nPhotos+i;
    ColIdx((nPhotos*nPixels)*2+i*3-2,1)=i;
    val((nPhotos*nPixels)*2+i*3-2,1)=1;
    
    RowIdx((nPhotos*nPixels)*2+i*3-1,1)=nPixels*nPhotos+i;
    ColIdx((nPhotos*nPixels)*2+i*3-1,1)=i+2;
    val((nPhotos*nPixels)*2+i*3-1,1)=1;
    
    RowIdx((nPhotos*nPixels)*2+i*3,1)=nPixels*nPhotos+i;
    ColIdx((nPhotos*nPixels)*2+i*3,1)=i+1;
    val((nPhotos*nPixels)*2+i*3,1)=-2;

   % coeffMtrx(nPixels*nPhotos+i,i-1)=1;
   % coeffMtrx(nPixels*nPhotos+i,i+1)=1;
   % coeffMtrx(nPixels*nPhotos+i,i)=-2;
end
coeffMtrx=sparse(RowIdx,ColIdx,val,nPixels*nPhotos+254,nPixels+256);
answer=coeffMtrx\lnT_gz;