function [betapulse,MSPE,NormModelSE,NaiveSE,RXbar,GXbar] = extract(imfile,imfil2,file,header)
%function [betapulse,MSPE,NormModelSE,NaiveSE,Flag,FlagGP] = extract(imfile,imfil2,file,header)
% EXTRACT Finds three measures for the relative signal for each spot
% K is the green channel and K2 is the red channel
% file is a file containing spot placement info, format 'blah.txt'
% header specifies the existence of column labels (1=TRUE, 0=FALSE)
% [betapulse,MSPE,NormModelSE,NaiveSE,Rxbar,Gxbar] = extract('2002-02-13_100-vs-50_replicate2_635_nm.tif','2002-02-13_100-vs-50_replicate2_532_nm.tif','2002-02-13_100-vs-50_rep2.txt',1);

red=imread(imfile);
green=imread(imfil2);
p=floor(19/2);
red=double(red); , green=double(green);
green(green==0)=1;

% For now we call this red out of laziness, even though its the log-ratio.
red = log2(red./green);
clear green

if header==1
[id,BlockID,Spotx,Spoty,flag,gridx,gridy,diamx,diamy,GXbar,GSE,GMedian,GBackest,GN,RXbar,RSE,RMedian,RBackest,RN,genName,ID,X,Y,Diam,F635Med,F635Mean,F635SD,B635Med,B635Mean,B635SD,F532Med,F532Mean,F532SD,B532Med,B532Mean,B532SD,B532_2SD,FPix,BPix,LogRatio,Flags]=textread(file,'%d%d%d%d%s%d%d%d%d%9.2f%9.2f%9.2f%9.2f%3d%9.2f%9.2f%9.2f%9.2f%3d%s%s%d%d%d%d%d%d%d%d%d%d%d%d%d%d%d%d%d%d%9.3f%d','delimiter','\t','headerlines',1);
names=textread(file,'%s',40,'delimiter','\t');
elseif header==0
[id,BlockID,Spotx,Spoty,flag,gridx,gridy,diamx,diamy,GXbar,GSE,GMedian,GBackest,GN,RXbar,RSE,RMedian,RBackest,RN,genName,ID,X,Y,Diam,F635Med,F635Mean,F635SD,B635Med,B635Mean,B635SD,F532Med,F532Mean,F532SD,B532Med,B532Mean,B532SD,B532_2SD,FPix,BPix,LogRatio,Flags]=textread(file,'%d%d%d%d%s%d%d%d%d%9.2f%9.2f%9.2f%9.2f%3d%9.2f%9.2f%9.2f%9.2f%3d%s%s%d%d%d%d%d%d%d%d%d%d%d%d%d%d%d%d%d%d%9.3f%d','delimiter','\t');
end

%clear GXbar GSE GMedian GBackest GN RXbar RSE RMedian RBackest RN genName ID X Y Diam F635Med F635Mean F635SD B635Med B635Mean B635SD F532Med F532Mean F532SD B532Med B532Mean B532SD B532_2SD FPix BPix LogRatio
clear GSE GMedian GBackest GN RSE RMedian RBackest RN genName ID X Y Diam F635Med F635Mean F635SD B635Med B635Mean B635SD F532Med F532Mean F532SD B532Med B532Mean B532SD B532_2SD FPix BPix LogRatio

nspots=length(gridx);
nnames=length(names);

iter=6608;
%iter=6144;

MSPE=zeros(1,iter);
betapulse=zeros(2,iter);
modelSE=zeros(1,iter);
NormModelSE=zeros(1,iter);
NaiveSE=zeros(1,iter);

Flag=char(flag); , poo=Flag(:);
Flag=poo(1:iter); , poo=Flag;
Flag=int8(poo);
% No flag and 'Artifact Noise'
Flag(poo==111)=0; , Flag(poo==65)=1;
% 'Irregular Ellipse' and 'Low Expression'
Flag(poo==73)=2; , Flag(poo==76)=3;
% 'No Detection' and 'Manual Spot'
Flag(poo==78)=4; , Flag(poo==77)=5;
Flag=double(Flag);

FlagGP = Flags;

%block=384;
for i=1:iter
i
    workred=red(gridy(i)-p:gridy(i)+p,gridx(i)-p:gridx(i)+p);

% Background calculation
%   workn=length(workred(:));
%        [numbin binval]=hist(workred(:),workn);
%        histpair=[numbin' binval'];
%        maxmode=histpair(histpair(:,1)==max(numbin),:);
%        bkr(i)=min(maxmode(:,2));
%   C=bkr(i)-min(workred(:));
%   workred=workred-bkr(i);

% Spot indicator
   [y0 x0] = size(workred);
   ax=diamx(i)/2;
   ay=diamy(i)/2;
   xx = (1:x0) - p;
   yy = (1:y0) - p;
   xx = (xx .* xx) / (ax * ax);
   yy = (yy .* yy) / (ay * ay);
   xx = repmat(xx,y0,1);
   yy = repmat(yy',1,x0);
   ind = xx + yy;
   inside = (ind <= 1);

mean(mean(workred(ind <= 1)))

   n=length(workred);   
   n=n(1);
   [beta, MSP] = spottest(n,workred,inside);

MSPE(i) = MSP;
betapulse(:,i) = beta;

% Parametric computation

result=spotse(workred,inside);
NaiveSE(i) = result(2);
NormModelSE(i) = result(3);
%NormModelSE(i) = 0;
end