function [newimage,dosv,dosi,dosj]=dosub(oldimage)
% [newimage,dosv,dosi,dosj]=DOSUB(oldimage)
%
% Takes a multispectral image (n x m x z) and turns it into a new
% multispectral image "hard-thresholded" to rendered the overall darkest
% object exactly zero, and all the other image pixes correspondingly, and
% suitably, adjusted. 
%
% INPUT:
%
% oldimage    The original n by m by z-channel image
%
% OUTPUT:
%
% newimage    The new n by m by z-channel image
% dosv        The z-channel darkest object
% dosi,dosj       The row/column number darkest object 
%
% EXAMPLE 1:
% 
% clear all; close all; figure(1)
% I=imread('pears.png');
% subplot(221);
% imshow(I)
% subplot(222);
% imshow(dosub(I))
% [nI,v,i,j]=dosub(I);
% subplot(223);
% imagesc(sum(double(I),3)); axis image; colorbar
% hold on; plot(j,i,'yo') ; hold off
% colormap(gray)
% subplot(224);
% imagesc(sum(double(nI),3)); axis image; colorbar
% hold on; plot(j,i,'yo') ; hold off
%
% EXAMPLE 2:
%
% clear all; close all
% load eb_enotre.mat
% ni=round(size(enotre_drone.ims.specdata,1)/3);
% mi=round(size(enotre_drone.ims.specdata,2)/3);
% ortho=enotre_drone.ims.specdata(ni:2*ni,mi:2*mi,:);
% oXX=enotre_drone.ims.XX(ni:2*ni,mi:2*mi,:);
% oYY=enotre_drone.ims.YY(ni:2*ni,mi:2*mi,:);
% figure(1)
% ah(1)=subplot(121);
% imagesc([oXX(1) oXX(end)],[oYY(1) oYY(end)],rapideya(ortho)); axis image xy
% ah(2)=subplot(122);
% imagesc([oXX(1) oXX(end)],[oYY(1) oYY(end)],rapideya(dosub(ortho))); axis image xy
% [nI,v,i,j]=dosub(ortho);
% hold on; plot(oXX(i,j),oYY(i,j),'yo') ; hold off
% 
% TESTED ON:
%
% 9.0.0.341360 (R2016a)
%

% Last modified by fjsimons-at-alum.mit.edu, 12/04/2019
% Find the darkest object, both where it is, and what it is
[doss,dosil]=min(sum(double(reshape(oldimage,prod(size(oldimage))/size(oldimage,3),[])),2));
% Return the row and column number...
[dosi,dosj]=ind2sub(size(oldimage(:,:,1)),dosil);
% ... and the image values of the darkest object
dosv=oldimage(dosi,dosj,:);
% Perform the thresholding. Note that this is only one of the many, many
% ways to do dehazing... and it's not subtracting the darkest object in
% every band separately, but finds the the darkest overall object and
% then subtracts the actual value at that position from every channel,
% whether that value itself was the minimum over the band or not. Note
% that we stay in the image variable class - uint16 bottoms out at zero
for index=1:size(oldimage,3)
  newimage(:,:,index)=oldimage(:,:,index)-oldimage(dosi,dosj,index);
end
  
