function result=spotse(y,indicator);
% SPOTSE assumes normality to compute a spot variance measure
% ans=spotse(y,indicator);
% y is a grid of pixel data
% indicator is spot membership

rho=0.5;
maxit=20;
epsilon=0.0001;

[rown,coln]=size(y);
signal=sum(sum(y.*indicator))/sum(sum(indicator))-sum(sum(y.*(1-indicator)))/sum(sum(1-indicator));
indicator=double(indicator);
meanvalue=indicator;n0=sum(sum(1-indicator));n1=sum(sum(indicator));
meanvalue(indicator==0)=sum(sum(y.*(1-indicator)))/n0;
meanvalue(indicator==1)=sum(sum(y.*indicator))/n1;
residuals=y-meanvalue;
res_vec=reshape(residuals,rown*coln,1);
sigma=std(res_vec);sigma2=sigma*sigma;

  sub_l=repmat((1:rown),1,rown)';
  sub_j=repmat((1:rown),rown,1); , sub_j=sub_j(:);
  sub_k=sub_l; , sub_m=sub_j;
  Sub_j=repmat(sub_j,1,rown^2);
  Sub_k=repmat(sub_k',rown^2,1);
  Sub_l=repmat(sub_l,1,rown^2);
  Sub_m=repmat(sub_m',rown^2,1);
%  dis=sqrt((Sub_j-Sub_l).^2+(Sub_k-Sub_m).^2);
  dis=max(abs(Sub_j-Sub_l),abs(Sub_k-Sub_m));
  temp=(y(sub_j,sub_l)-meanvalue(sub_j,sub_l)).*(y(sub_k,sub_m)-meanvalue(sub_k,sub_m))/sigma2;
% Eliminating correlation after distances between pairs of one
% The temp variable must be zeroed out first
temp(dis>1)=0;
dis(dis>1)=0;
  dis=dis(:); , temp=temp(:);
clear Sub_*  

% This takes the longest
converg=0;it=1;
while converg==0,
%    ls = sum((temp-rho.^dis).^2);
    du = sum((temp-rho.^dis).*(rho.^(dis-1)).*dis);
    ddu = sum((temp.*(dis-1).*(rho.^(dis-2)).*((dis-2)>0)-(2*dis-1).*(rho.^(2*dis-2))).*dis);
    rho1=rho-du/ddu;  
if rho1 < 0
warning('Correlation estimate is negative')
end
%    disp([it rho rho1]);    
    converg=abs(rho1-rho)<epsilon;
    rho=rho1;
    it=it+1;
    if it>maxit,break;end;
end;

c_matrix=indicator*n0-(1-indicator)*n1;
tmpC = c_matrix(sub_j,sub_l).*c_matrix(sub_k,sub_m);
clear sub_*

tmpC = tmpC(:);
tempC = sum(tmpC(dis>0).*(rho.^dis(dis>0)));

SE_est=1/(n0*n1) * sigma * sqrt(sum(sum(c_matrix.*c_matrix))+tempC);
Naive_SE=sigma*sqrt(1/n0+1/n1);
% disp(['Estimated means=' num2str(signal) ';  Estimated Standard Error=' num2str(SE_est) ';  Naive Estimate=' num2str(Naive_SE)]);
result=[signal Naive_SE SE_est];