clear all;
   
  %starvärden för vissa variabler
  count=512*512;
  bild=0;bild2=0;bild3=zeros(512);
  bildtot=0;bildtot2=0;bildtot3=0;
  h(1:2000)=0;hh(1:2000)=0;hhh(1:2000)=0;hhhh(1:2000)=0;
  m=0;mm=0;mmm=0;
  c=0;
 
  % mallar för utsmetade fotoner
  mat(:,:,1)=[1 0 1 1;0 1 0 1;1 0 1 1;1 1 1 1];mat(:,:,2)=[1 0 1 1;0 1 0 1;0 1 0 1;1 0 1 1];mat(:,:,3)=[1 0 0 1 ;0 1 1 0 ;1 0 0 1;1 1 1 1];mat(:,:,4)=[1 0 1 1 ;0 1 0 1;0 1 1 0;1 0 0 1];
  mat(:,:,5)=[1 0 0 1;0 1 1 0;0 1 0 1 ;1 0 1 1];mat(:,:,6)=[1 0 0 1;0 1 1 0;1 0 1 0;1 1 0 1];mat(:,:,7)=[1 0 0 1; 0 1 1 0; 0 1 1 0;1 0 0 1];mat(:,:,8)=[1 1 0 1;1 0 1 0;0 1 1 0;1 0 0 1];
  
  
  % inläsning av bilden
  f=fopen('insb1x60.spe');
  fread(f, 4100, 'uint8'); % slänger skräpet i början
  
 
  %whilesats, behövlig om man har flera bilder i samma fil
  while count==512*512 
    	  [i, count]=fread(f, [512 512], 'uint16'); 
   
      if count==512*512 
          bild=i;
          clear a foton n x;
          a=max(max(bild));
          
          if a>m % m bestämmer längden på slutgiltiga plotten
              m=a;
          end
          
          % Gör ett histogram över obehandlad delbild
          x=0:3:a;
          n=hist(bild,x);
          foton=sum(n');
          a=round(a/3);
          h(1:a-1)=h(1:a-1)+foton(1:a-1); % lägger ihop delhistogrammen till ett enda
          
          bildtot=bildtot+bild; % lägger ihop delbilderna till en enda
         
                    
          %Gör en provisorisk bakgrundsberäkning:
          bak=0;  
          a=(max(bild));
          a=max(a);
              
           x=0:3:a;
           n=hist(bild,x);
           foton=sum(n');
           bak=find(foton==max(foton));  %letar upp bakgrundstoppen                 
         
           %Titta på det som ej verkar vara bakgrund (>bakgrunden+6)
           mall=bild>(bak*3+6); 
           bild2=mall.*bild;
           [u,v,varde]=find(bild2); % (u,v) är koordinater och varde är värdet
          
           % Gå igenom alla pixlar som är fotoner 
           % Detta tar någon minut...
         
           bild3(1:512,1:512)=0;
           for k=1:(length(u))
               uu=u(k);vv=v(k);
               
               if uu==1 | vv==1 | uu==512 | vv==512
                   bild2(uu,vv)=(bak*3+6);% om det är en pixel i kanten så kastas den (blir bakgrund)
                   %lokala bakgrunden får ett approximativt värde och p sätts till noll
                   medel=bak*3+6;
                   p=0;
               else
                              
              %Titta på åtta pixlar runt omkring fotonen     
              bildsort=(bild(uu-1:uu+1,vv-1:vv+1));
              %lägg i vektor och sortera
               bildvek=[bildsort(1,:) bildsort(2,:) bildsort(3,:)];
               s=sort(bildvek);
               %ta medelvärdet av de lägsta (bör bli lokal bakgrund)
               p=find(s<bak*3+6);
               p=length(p);
               end %if-end
               
               % om området har 'bakgrundspixlar' så beräknas lokal bakgrund
               if p>0 
               medel=round(mean(s(1:p)));
               end % if-end
               
               %bakgrunden tas bort, beräknad där det varit möjligt och den approximativa annars
               bild2(uu,vv)=bild2(uu,vv)-medel;   
                 
           end % for-end
          
           
           a=max(max(bild2));
           if a>mm %mm bestämmer längden på slutgiltiga plotten
              mm=a;
           end
          
           %gör histogram efter bakgrundsubtraktion
           xx=0:3:a;
           nn=hist(bild2,xx);
           foton2=sum(nn');
           a=round(a/3);
           hh(1:a-1)=hh(1:a-1)+foton2(1:a-1);
          
           bildtot2=bildtot2+bild2;
         
          % lägg ihop fotondelar som spritts ut
          % gå igenom de pixlar som är fotoner (större än bakgrund)

            for j=1:length(u);
                uu=u(j);vv=v(j);
                
                 
                if uu>1 & uu<511 & vv>2 & vv<511 % fotoner i kanten får kastas... kolonn måste vara större är 2 för fallet mat8
                if bild2(uu,vv-1)==0 & bild2(uu-1,vv)==0 % så inte fotonen hanteras mer än en gång
                   
                    pp=1; % startvärde för while-loopen
                   
                   while pp<=8
                       if pp<8
                           if bild2(uu-1:uu+2,vv-1:vv+2).*mat(:,:,pp)==bild2(uu-1:uu+2,vv-1:vv+2)
                               break
                           end
                       else
                           if bild2(uu-1:uu+2,vv-2:vv+1).*mat(:,:,pp)==bild2(uu-1:uu+2,vv-1:vv+2)
                               break
                           end
                       end
                       pp=pp+1;
                   end
                 
                   if pp>1 & pp<8
                       bilds=sum(bild2(uu:uu+1,vv:vv+1));
                       bilds=sum(bilds');
                       bild3(uu,vv)=bilds;
                   
                   elseif pp==8
                       bilds=sum(sum(bild2(uu:uu+1,vv-1:vv)));
                       bild3(uu,vv)=bilds;
                                     
                  elseif pp==1
                   bild3(uu,vv)=bild2(uu,vv);
                  
                  end
                           
               end%if-end
               end%if-end
           end%for-end 
               
              a=max(max(bild3));
              if a>mmm 
              mmm=a;
              end

              %Gör histogram över modifierad delbild
              xxx=0:3:a;
              nnn=hist(bild3,xxx);
              foton3=sum(nnn');
              a=round(a/3);
              hhh(1:a-1)=hhh(1:a-1)+foton3(1:a-1);
          
         bildtot3=bildtot3+bild3;
         
               c=c+1 % visar hur många bilder som behandlats 
             
         end % if-end
  end % while-end
  
  
  %Plotta de tre olika totalhistogrammen
  
  figure(1);
  x=0:3:m;
  plot(x(1:round(m/3)),h(1:round(m/3)));
  
  figure(2);
  xx=0:3:mm;
  plot(xx(10:round(mm/3))*0.020038,hh(10:round(mm/3)));
  
  figure(3);
  xxx=0:3:mmm;
  plot(xxx(10:round(mmm/3))*0.020038,hhh(10:round(mmm/3)));
  
