function snr_curve(raw1, raw2, black_level, white_level, fit_model)
    % extracts the SNR curve from two test images, taken with identical settings
    % from the same static scene; method inspired from (but not identical to)
    % http://www.clarkvision.com/articles/evaluation-1d2/ 
    %
    % parameters:
    % - raw1, raw2: CR2 file names, to be read via dcraw
    %               (so, make sure dcraw is in your PATH)
    % - black_level, white_level: find them with exiftool
    %              (white level may be also found with prctile(im(:),99.99))
    % - fit_model: optionally fit a noise model, estimating how much
    %              is read noise and how much is Poisson shot noise,
    %              and figuring out the gain (e/DN) and the full-well capacity
    % 
    % (c) 2014 a1ex. GPL.

    disp('reading images...')
    im1 = read_cr2(raw1, 0);
    im2 = read_cr2(raw2, 0);
    
    imdif = im1 - im2;
    imavg = (im1 + im2) / 2;
    
    signal = [];
    noise = [];
    
    [H, W] = size(im1);
    disp('sampling...')
    
    for k = 1:10000
        x = round(50 + rand() * (W - 100));
        y = round(50 + rand() * (H - 100));
        [avg, unused] = sample_patch(imavg, x, y, 8);
        if avg < white_level - 300
            [unused, stdev] = sample_patch(imdif, x, y, 8);
            stdev = stdev / sqrt(2);
            signal(k) = log2(max(1, avg - black_level));
            noise(k) = log2(stdev);
        end
    end
    
    if fit_model
        disp('fitting...')
        model = fminsearch(@(x) check_noise_model(x, signal, noise), [8 1]);
        read_noise = model(1);
        gain = model(2);
        
        read_noise
        gain
        
        model_signal = linspace(0, log2(white_level - black_level), 100);
        dn = 2.^model_signal;
        electrons = dn * gain;
        shot_snr = sqrt(electrons);
        shot_noise = dn ./ shot_snr;
        model_noise = log2(sqrt(read_noise.^2 + shot_noise.^2));
        full_well = round((white_level - black_level) * gain);
        full_well
        
        plot(signal, signal-noise, '.b', model_signal, model_signal-model_noise, '-r', 'linewidth', 3);
        legend('Measured SNR', 'Fitted SNR model', 'location', 'southeast')
    else
        plot(signal, signal-noise, '.b');
        legend('Measured SNR', 'location', 'southeast')
    end
    grid on
    max_signal = log2(white_level - black_level);
    axis([max_signal-15 max_signal -2 10])
    xlabel('Input signal (EV), showing 15 stops below white level');
    ylabel('SNR, from -2 EV to +10 EV');
end

function err = check_noise_model(model, input_signal, measured_noise)
    read_noise = model(1);
    gain = model(2);
    
    dn = 2 .^ input_signal;
    electrons = dn * gain;
    shot_snr = sqrt(electrons);
    shot_noise = dn ./ shot_snr;
    combined_noise = sqrt(read_noise.^2 + shot_noise.^2);
    delta = measured_noise - log2(combined_noise);
    err = norm(delta);
end

function [avg, stdev] = sample_patch(im, x, y, radius)
    % sample channels of the same color from a Bayer grid
    patch = im(y-radius : 2 : y+radius, x-radius : 2 :x+radius);
    avg = mean(patch(:));
    stdev = std(patch(:));
end

function im = read_cr2(filename, also_ob)
    if also_ob
        dcraw = "dcraw -c -4 -E";
    else
        dcraw = "dcraw -c -4 -D";
    end
    system(sprintf('%s %s > tmp.pgm', dcraw, filename));
    im = double(imread('tmp.pgm'));
    system('rm tmp.pgm');
end