function camData= DXOdataFit()

%fit curve to DXOmark sensor SNR data to estimate

%pixel read noise, saturation count, and standard deviation of gain error

%Noise model is: noise^2 = n^2 + m*x + G * m^2 * x^2

%x : fraction of saturation count

%G : pixel gain variance

%m : saturation count

%n : read noise

% Outputs:

% camData.name : Camera data file name

% camData.ISO : Manufacturers ISO values

% camData.readNoise : pixel read noise in electrons (n)

% camData.eSat : saturation count in electrons (m)

% camData.PRNU : standard deviation of noise gain error (sqrt(G)

% camData.res % residual (difference) of fit and data

plotResidues= 0; % plot difference between fit and data?;

fprintf(1,'\n* Sensor Paramaters Derived From DXOmark.com Data *');

fprintf(1,'\n\nISO: Manufacturers ISO rating');

fprintf(1,'\nm: Pixel saturation count in electrons');

fprintf(1,'\nn: Pixel read noise in electrons');

fprintf(1,'\ng%%: Standard deviation of gain error in percent');

fprintf(1,'\nres: Worst case difference between fitted SNR curve in dB and SNR data points in dB');

fprintf(1,'\n');

[xISO,SN,ISO,cam]= dataread;

iCam= 1;

while ~isempty(ISO)

fprintf(1,['\n' cam]);

fprintf(1,'\nISO:');

fprintf(1,' %5.0f',ISO);

if plotResidues

figure(1);

clf;

CO= get(gca,'ColorOrder');

end;

Niso= length(ISO);

n= zeros(1,Niso);

G= zeros(1,Niso);

m= zeros(1,Niso);

res= zeros(1,Niso);

for i1= 1:Niso

% for balanced weighting of the points being fitted

% fit uses x*noise^2/(m*x)2 = n^2/m^2 * 1/x + 1/m + G * x

xp= 1e-2*xISO{i1}; % x values to plot

x= xp(2:end-1); % x values to fit

SNRp= SN{i1};

SNR= SNRp(2:end-1);

N2SR= x .* 10.^(-SNR./10); % x * noise^2 / (m*x)^2

Ap= [xp.^-1 xp.^0 xp.^1];

% Fit the data

A= [x.^-1 x.^0 x.^1];

c= (A' *A)\(A' *N2SR); % least squares solution of c = A\N2SR

G(i1)= c(3);

n(i1)= sqrt(c(1))./c(2);

m(i1)= 1/c(2);

SNRfit= -10*log10(Ap*c./xp);

res(i1)= max(abs(SNRfit(2:end-1)-SNR));

if plotResidues

semilogx(xp(2:end),SNRfit(2:end)-SNRp(2:end),'Color',CO(mod(i1-1,size(CO,1))+1,:));

hold on

end

end

if plotResidues

hold off

xlabel('Signal');

ylabel('Difference of dB');

title([cam ' Fit Residuals']);

legend(num2str(ISO'),2);

drawnow;

end

camData(iCam).name= cam;

camData(iCam).ISO= ISO;

camData(iCam).eSat= m;

camData(iCam).readNoise= n;

camData(iCam).PRNU= sign(G).*sqrt(abs(G));

camData(iCam).res= res;

fprintf(1,'\n m = ');

fprintf(1,'%5.0f ',m);

fprintf(1,'\n n = ');

fprintf(1,'%5.2f ',n);

fprintf(1,'\ng%% = ');

fprintf(1,'%5.2f ',sign(G).*sqrt(abs(G))*100);

fprintf(1,'\nres= ');

fprintf(1,'%5.3f ',res);

fprintf(1,'\n');

[xISO,SN,ISO,cam]= dataread;

iCam= iCam+1;

end

end

function [xISO,SN,ISO,cam]= dataread()

[fname,pname]= uigetfile('*.php','select snrdetail source file');

xISO= [];

SN= [];

ISO= [];

cam= [];

if ~(fname==0)

[~,cam]= fileparts(fname);

fid=fopen(fullfile(pname,fname),'r');

done= 0;

i2= 1;

while ~done;

a{1}=[];

while(~done && isempty(a{1}))

C = textscan(fid, '%s',1,'delimiter','>');

a=(strfind(C{1}, 'dataSet'));

if isempty(a)

done= 1;

end

end;

if ~done

beginPt= ftell(fid);

a=(strfind(C{1}, 'ISO'));

str= C{1}{1};

ISO= [ISO, sscanf(str(a{1}+4:end),'%f')];

a{1}=[];

i1=0;

end;

while(~done && isempty(a{1}))

C = textscan(fid, '%s',1,'delimiter','>');

a=(strfind(C{1}, 'set'));

if isempty(a)

done= 1;

else

if ~isempty(a{1})

i1= i1+1;

end

a=(strfind(C{1}, '/dataSet'));

end

end;

if ~done

fseek(fid,beginPt,'bof');

C = textscan(fid, '%s',5*i1,'delimiter','''');

C= C{1};

C= reshape(C,[5 i1]);

y= str2num(char(C{2,:}));

x= str2num(char(C{4,:}));

% Old code to remove bad points from old style files disabled %%%%%%%%%

% xt= [.01:.01:.09 .1:.1:.5 .7 1:1:9 10:1:14 15:5:35 42:1:255];

% x(2:end)= xt(end-length(x)+2:end); % fix the bad x values

% xISO{i2}= x*100/255;

% SN{i2}= y;

xISO{i2}= x(2:end)*100/255;

SN{i2}= y(2:end);

i2= i2+1;

a{1}=[];

while(~done && isempty(a{1}))

C = textscan(fid, '%s',1,'delimiter','>');

if isempty(a)

done= 1;

else

a=(strfind(C{1}, '/dataSet'));

end

end

end

end

fclose(fid);

end

end