applyRawFlatFrame.m

%%
%% Applies a raw flat frame to a raw image using zero blur, which allows the
%% removal of sensor dust spots (and also vignetting correction) at the expense
%% of introducing some noise from the flat frame
%%
%% Parameters:
%%
%% imgFilename        - Filename of image to correct, generated via LibRaw's unprocessed_raw.exe
%% flatFieldFilename  - Filename of flat-field image, generated via LibRaw's unprocessed_raw.exe
%% outputFilename     - Filename to hold output bayer PGM, which caller needs to convert toc
%%                      raw via pgm2dng.exe
%% imgDimensions      - Dimensions of image to use. Example: [4024 6018] for 6018x4024
%%
%% Output:
%%
%% PGM bayer-interlaced file, corrected for sensor dust via flat-frame
%%
function applyRawFlatFrame(imgFilename, flatFieldFilename, outputFilename, imgDimensions)

  cfaPatternStr = "RGGB";

  function imgData = loadImgData(filename)
    imgData = double(imread(filename));
    imgData = resize(imgData, imgDimensions);
  end

  function avg = calcChannelCenterAverage(channelData)
    avg = mean(channelData(floor(end/2-128):floor(end/2+128), floor(end/2-128):floor(end/2+128))(:));
  end

  %% read target and flat-field raw data generated by LibRaw's unprocessed_raw.exe
  imgData = loadImgData(imgFilename, imgDimensions);
  imgDataFlatField = loadImgData(flatFieldFilename, imgDimensions);

  %% split the raw data into its separate RGGB channels
  [r, g1, g2, b] = rawBayerToChannels(cfaPatternStr, imgData);
  [rFlat, g1Flat, g2Flat, bFlat] = rawBayerToChannels(cfaPatternStr, imgDataFlatField);

  %% generate flat-field divisor tables for each color channel
  rDivisors  = rFlat  ./ calcChannelCenterAverage(r);
  g1Divisors = g1Flat ./ calcChannelCenterAverage(g1);
  g2Divisors = g2Flat ./ calcChannelCenterAverage(g2);
  bDivisors  = bFlat  ./ calcChannelCenterAverage(b);

  %% apply flat-field divisors to image data
  r  = r  ./ rDivisors;
  g1 = g1 ./ g1Divisors;
  g2 = g2 ./ g2Divisors;
  b  = b  ./ bDivisors;

  %% put the modified channels back into a bayer pattern
  imgDataOut = channelsToRawBayer(cfaPatternStr, r, g1, g2, b);

  %% write a PGM with the new data
  imwrite(uint16(imgDataOut), outputFilename);
end