data = double(imread('300D_res.ppm'));

%for simplicity, just use green channel
image = data(:,:,2);

MTF_300d = zeros(13,1);
frequencies = [1 2 3 4 5 6 8 10 12 14 16 18 20];


%1 through 6
curr_region = image(677, 1190:1350);
MTF_300d(1) = (max(curr_region)-min(curr_region))/(max(curr_region)+min(curr_region));
curr_region = image(800, 1230:1312);
MTF_300d(2) = (max(curr_region)-min(curr_region))/(max(curr_region)+min(curr_region));
curr_region = image(926, 1240:1300);
MTF_300d(3) = (max(curr_region)-min(curr_region))/(max(curr_region)+min(curr_region));
curr_region = image(1052, 1250:1290);
MTF_300d(4) = (max(curr_region)-min(curr_region))/(max(curr_region)+min(curr_region));
curr_region = image(1178, 1255:1287);
MTF_300d(5) = (max(curr_region)-min(curr_region))/(max(curr_region)+min(curr_region));
curr_region = image(1301, 1257:1286);
MTF_300d(6) = (max(curr_region)-min(curr_region))/(max(curr_region)+min(curr_region));

%8 through 20
curr_region = image(802, 1780:1820);
MTF_300d(7) = (max(curr_region)-min(curr_region))/(max(curr_region)+min(curr_region));
curr_region = image(885, 1783:1817);
MTF_300d(8) = (max(curr_region)-min(curr_region))/(max(curr_region)+min(curr_region));
curr_region = image(967, 1786:1813);
MTF_300d(9) = (max(curr_region)-min(curr_region))/(max(curr_region)+min(curr_region));
curr_region = image(1050, 1788:1813);
MTF_300d(10) = (max(curr_region)-min(curr_region))/(max(curr_region)+min(curr_region));
curr_region = image(1132, 1790:1811);
MTF_300d(11) = (max(curr_region)-min(curr_region))/(max(curr_region)+min(curr_region));
curr_region = image(1215, 1791:1809);
MTF_300d(12) = (max(curr_region)-min(curr_region))/(max(curr_region)+min(curr_region));
curr_region = image(1294, 1792:1808);
MTF_300d(13) = (max(curr_region)-min(curr_region))/(max(curr_region)+min(curr_region));