flatinfo = dir( fullfile(workdir, 'flat*') );
flatlist = fullfile( workdir, {flatinfo.name} );
nflat = numel(flatlist);
darkinfo = dir( fullfile(workdir, 'dark*') );
darklist = fullfile( workdir, {darkinfo.name} );
ndark = numel(darklist);
flat = zeros(ysize,xsize);
for K = 1 : nflat
imtemp = imresize( imread(flatlist{k}), [ysize, xsize] );
flat = flat + double(imtemp);
end
flat = flat ./ nflat;
dark = zeros(ysize,xsize);
for K = 1 : ndark
imtemp = imresize( imread(darklist{k}), [ysize, xsize] );
dark = dark + double(imtemp);
end
dark = dark ./ ndark;
flatsub = flat - dark;
flatsub(flatsub <= 0) = 1;
But are you sure that you want to set the zero / negative areas to 1 and not to 0? To me it would make more sense to use
flatsub = max(0, flat - dark);
which would use 0 for places that would have been negative.