# A very crude tool for applying non-uniform gain to integrating a series of image data
# Makes a normal cubic spline based on user-supplied control points to model the gain function for each region
# Applies the average gain for the frame interval sampled from that integration function 
# Selects gain functions based on a set of regions specfied as an input-sized pgm with one gray per region/function.

#### This part is setup ####
#Image Input Setup
fr=(1/240); #framerate of input sequence, in seconds (not sensed, video package broken on 5.x)
path="20190930_205109.mp4"; #Path to source file

#Read in an appropriate user-created mask
mask=imread('PenguinRockMask.pgm');
mask=cast(mask,'double');
#Mask values defined in the mask
# would be better to get from mask automatically, 
# but I keep seeing border pixels w/ garbage
gray=[255,128];
fns=length(unique(gray));

#Integration Functions Setup 
pts= [100, 200]; 
x=0:(max(pts)-1);
y=zeros(fns,max(pts));
y(1,1:pts(1))=y=[0 0 0 0 0 0 0 0 0 10 30 50 30 10 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 10 30 50 30 10 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 ];
y(2,1:pts(2))=[0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 400 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0  0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0  0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0  0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0];

#### End of configuration ####

#Crude changem cribbed from, 
# https://stackoverflow.com/questions/11952037/replace-values-in-matrix-with-other-values
# because I'm not using a commercial toolbox for this
function mapout = changem(Z, newcode, oldcode)
% Identical to the Mapping Toolbox's changem
% Note the weird order: newcode, oldcode. I left it unchanged from Matlab.
    if numel(newcode) ~= numel(oldcode)
        error('newcode and oldcode must be equal length');
    end

    mapout = Z;

    for ii = 1:numel(oldcode)
        mapout(Z == oldcode(ii)) = newcode(ii);
    end
end

#Show the integration functions (comment to hide)
for dfs=1:fns
  figure(dfs);
  xspline=0:0.01:pts(dfs); 
  yspline=spline(x(1:pts(dfs)),y(dfs,1:pts(dfs)),xspline);
  ylabel("gain");
  xlabel("fractional time");
  plot(xspline,yspline,"g-",x(1:pts(dfs)),y(dfs,1:pts(dfs)),"b+");
endfor


#Have ffmpeg process video for frames
# In a better world, I'd use the video package, but it's currently broken on 5.x
# https://savannah.gnu.org/bugs/?51344
# That would give automatic frame-rate extraction, but nooo
##ffmpegline=["ffmpeg -i ", path ," -vcodec ppm img%03d.ppm"]
##system(ffmpegline);

t=0; #time
fc=0; #number of frames ingested, maintained seperately because I'm lazy
sumimg=0;

#for each input image
files = dir('img*.ppm');
numimg = length(files);
texp = fr*numimg; 
for file = files'
  #read in image
  img=imread(file.name);
  #Promote for range
  img= im2double(img);
  #compute gains for that frame 
  for gfs=1:fns
    #Still don't quite touch the last control point
    gain(gfs)=( ...
          interp1(x(1:pts(gfs)),y(gfs,1:pts(gfs)),(t/texp)*(pts(gfs)-1.01),"spline") + ...
          interp1(x(1:pts(gfs)),y(gfs,1:pts(gfs)),((t+fr)/texp)*(pts(gfs)-1.01),"spline"))/2;
    #printf("Frame %d gain %d\n", fc, gain);
  endfor
  #apply to frame
  
  
  gainmask=changem(mask,gain,gray);
  imgmod=img.*gainmask;
  #add frame to sum
  sumimg=imgmod+sumimg;
  fc=fc+1; #incriment frame count
  t=t+fr; #incriment time
endfor


#Divide by number of frames
finalimg=sumimg/fc;
#cleanup
##delete img*.ppm;
#save and display rendered image
imwrite(finalimg,"result.ppm");
figure(gfs+1);
imshow(finalimg);


