function [ output_args ] = fret_rescaling_reslicing_colormap_functions( turquoise, venus, draq5, NumberImages, scale, frequencyX, frequencyY, frequencyZ, resliced_directory, smoothing, iterations, steps, reduc, threshold)

%%=====================Turquoise and Venus Smoothing========================
%%Smooth Turquoise and Venus with "Smooth3" function
    turquoise = double(turquoise); %convert to double
    venus = double(venus); %convert to double

%% smoothing algorithems

%gaussian
if smoothing == 1
turquoise= smooth3(turquoise,'gaussian'); %takes doubles
venus = smooth3(venus,'gaussian'); %takes doubles
end

% 2d slice
if smoothing == 2
turquoise = aniso_2d_slices(turquoise,iterations, steps, reduc, 'Turquoise '); 
venus = aniso_2d_slices(venus,iterations, steps, reduc, 'Venus '); 
end

% 3d long
if smoothing == 3
turquoise = anisodiff3d(turquoise, iterations, steps, reduc, 2, [1 1 1]); 
venus = anisodiff3d(venus, iterations, steps, reduc, 2, [1 1 1]); 
end

% 3d fast
if smoothing == 4
turquoise = fast_anisodiff3d(turquoise, iterations, steps, reduc, 2, [1 1 1]); 
venus = fast_anisodiff3d(venus, iterations, steps, reduc, 2, [1 1 1]); 
end

    turquoise = uint16(turquoise); %turns to 16 bit class
    venus = uint16(venus);  %turns to 16 bit class
    
%%
%=====================Find Threshold==============================
%         
%=====================Turquoise and Venus Stretching=======================

    function [ ImageNEW ] = stretch(Image,scaling) %"Image is just turquoise or venus stack
        Image = permute(Image,[1 3 2]); %132,213,231,312,321
        
        num_images = length(Image(1, 1, :));
        [M,N] = size(Image(:, :, 1));
        ImageNEW = imrotate(zeros(M,ceil(N.*scaling),num_images),-90);
        for K=1:num_images
            %Slicing on X plane; must rotate if sliced on x plane
            ImagePass = imrotate(Image(:, :, K),-90); %rotate image 90 degrees clockwise
            %stretch the image
            [x,y] = size(ImagePass);
            toScale = x.*scaling;
            ImageNEW(:,:,K) = imresize(ImagePass,[toScale y]);                
            
        end
        ImageNEW = permute(ImageNEW,[3 2 1]);
        ImageNEW = imrotate(ImageNEW,90);
    end


turquoise = stretch(turquoise,scale);
venus = stretch(venus,scale);
draq5 = stretch(draq5,scale);

%=================FRET Calculations and Masking============================

%Convert images to double

        turquoiseAsDouble=double(turquoise);
        venusAsDouble=double(venus);
        draq5AsDouble=double(draq5);
      
    [mImage,nImage,oImage] = size(turquoise);
    FRETOut = zeros(mImage,nImage,oImage,'uint16');
    
     %Make mask from turquoise + venus and apply it to the FRETimage
        
        thresholdTV = turquoise + venus;

%This is for dilation/erosion of the following masks
        seClose = 5;
        seClose = strel('disk',seClose,0);
        
        maxthresholdTV = max(thresholdTV,[],3); %find max Thershold
        maxthresholdTV = multithresh(maxthresholdTV)*threshold;
        
    Waitbar_Handle = waitbar(0,'Please wait...','Name','Waiting...');
    for K=1:oImage
        
            waitbar(K/oImage,Waitbar_Handle,sprintf(strcat('"FRETing "',num2str(K/oImage*100),'%'),K,oImage))
        
        %Venus and Turquoise for a slice
          venus1=venus(:,:,K);
          turquoise1=turquoise(:,:,K);
    
            
        %thresholdTV = multithresh(thresholdTVTemplate);
        maskTV1 = thresholdTV(:,:,K) > maxthresholdTV; %make the mask
        maskTV1 = imfill(maskTV1,'holes'); %fill in the holes
        maskTV1 = imclose(maskTV1,seClose); %dilate and erode
        maskTV1 = imfill(maskTV1,'holes'); %fill in the holes again
        %FRET(~finalMaskTV_2) = 0; %apply the mask to FRET image
        
        %Make mask for DRAQ5 in order to remove it from FRET image
        midSlice = floor((oImage)./2);
        thresholddraq5 = draq5(:,:,midSlice);
        thresholddraq5 = multithresh(thresholddraq5);
        thresholddraq5 = draq5AsDouble(:,:,K) > thresholddraq5; %make the mask
        thresholddraq5 = imfill(thresholddraq5,'holes'); %fill in the holes
        thresholddraq5 = bwareaopen(thresholddraq5, 300);%remove small objects from DRAQ5 (Use to 1500 initially)
        thresholddraq5 = imclose(thresholddraq5,seClose); %dilate and erode
        thresholddraq5 = imfill(thresholddraq5,'holes'); %fill in the holes again
        %FRET(finalMaskDraq5_2) = 0; %apply the mask to FRET image
        
        %Mask total
        finalMaskdraq_inv = ~thresholddraq5;
        binaryMaskTotalFinal = maskTV1.*finalMaskdraq_inv; 
        
        %remove objects smaller than specified size from finalmask
        binaryMaskTotalFinal = ~binaryMaskTotalFinal;
        binaryMaskTotalFinal = bwareaopen(binaryMaskTotalFinal, 500);%remove small holes 
        binaryMaskTotalFinal = ~binaryMaskTotalFinal;
        binaryMaskTotalFinal = bwareaopen(binaryMaskTotalFinal, 500);%remove small objects

       
        %FRET equation is taken from Vigel et.al., 
        % 17.128 is correction factor calculated using Turquoise and Venus extinction coeffiencts @ 405 excitation
        
        FRET =   binaryMaskTotalFinal.*((venus1-turquoise1)./(turquoise1.*17.128+venus1).*(2.^16));

        %Convert back to uint16
        FRETOut(:,:,K) = uint16(FRET);
       
            %save masks
       outputFileName1 = strcat(resliced_directory,'\MASK_TV.tif');
        imwrite(maskTV1,outputFileName1, 'WriteMode', 'append',  'Compression','none');
       
       outputFileName1 = strcat(resliced_directory,'\MASK_DRAQ5.tif');
       imwrite(thresholddraq5,outputFileName1, 'WriteMode', 'append',  'Compression','none');
        
       outputFileName1 = strcat(resliced_directory,'\MASK_TOTAL_FINAL.tif');
       imwrite(binaryMaskTotalFinal,outputFileName1, 'WriteMode', 'append',  'Compression','none');
        
    end

    close(Waitbar_Handle)
    

 %%       
    
%=========================Reslicing and Saving=============================

    function [ satCount] = reslicing(Image,reslicing_directory,plane,frequency,name)
        new_directory=strcat(reslicing_directory,'\',num2str(frequency),'px','Resliced_',plane,name);
        mkdir(new_directory);
        
        ImageNEW = [];
        
        if 1 == strcmp('X',plane) %strcmp returns 1 for true, 0 for false
            Image = imrotate(permute(Image,[3 1 2]),90); %132-good (213,231, rotate 321 bottom-top)
            num_images = ceil(length(Image(1, 1, :))./frequency);
            [M,N] = size(Image(:, :, 1));
            ImageNEW = imrotate(zeros(M,N,num_images),-90);
            
        end
        
        if 1 == strcmp('Y',plane) %strcmp returns 1 for true, 0 for false
            Image = imrotate(permute(Image,[2 3 1]),-90); %321-good (213,)
            num_images = ceil(length(Image(1, 1, :))./frequency);
            [M,N] = size(Image(:, :, 1));
            ImageNEW = zeros(M,N,num_images);
        end
        
        if 1 == strcmp('Z',plane)
            num_images = ceil(length(Image(1, 1, :))./frequency);
            [M,N] = size(Image(:, :, 1));
            ImageNEW = zeros(M,N,num_images);
          
            
        end
        depth=length(Image(1, 1, :));
        satCount = [];%declare array for recording over and under saturation count for each image
        
        Waitbar_Handle = waitbar(0,'Please wait...','Name','Waiting...');
        for K=1:frequency:depth
            
            waitbar(K/depth,Waitbar_Handle,sprintf(strcat('Reslicing'),K,depth))
            
            %must rotate if sliced on x plane
            if 1 == strcmp('X',plane) %strcmp returns 1 for true, 0 for false
                ImageNEW(:, :, K) = imrotate(Image(:, :, K),-90); %rotate image 90 degrees clockwise
             waitbar(K/depth,Waitbar_Handle,sprintf(strcat('Reslicing X'),K,depth))
            end
            
            if 1 == strcmp('Y',plane)
               ImageNEW = Image;
             waitbar(K/depth,Waitbar_Handle,sprintf(strcat('Reslicing Y'),K,depth))
            end
            
            if 1 == strcmp('Z',plane)
               ImageNEW = Image;
   
            
            end
            
      
            [M,N] = size(ImageNEW(:,:,1));
            maskFRET = double(ImageNEW(:, :, K))./2.^16;
            addedValues = 0.0;
            n = 0; %count the number of added values
            for j = 1:M
                for k = 1:N
                    if maskFRET(j,k)>0
                        addedValues = addedValues + maskFRET(j,k);
                        n = n + 1;
                    end
                end
            end

            addedValues./n;
                     
                
            if strcmp(name,'FRET') %if FRET image, apply colormap before saving
                [FRETColor,overSat,overSatFrac,underSat,undrSatFrac] = grs2rgb(ImageNEW(:,:,K),0.60,0); %third parameter is R or A for Relative or Absolute

                imwrite(FRETColor,strcat(new_directory,'\',num2str(K),'.tif'),'tif');
                
                imwrite(uint16(ImageNEW(:,:,K)),strcat(new_directory,'\RawFRET_',num2str(K),'.tif'),'tif');
                
            else %else save as normal
                imwrite(ImageNEW(:, :, K),strcat(new_directory,'\',num2str(K),'.tif'),'tif');
            end
                       
        end
        

    end

reslicing(FRETOut,resliced_directory,'X',frequencyX,'FRET');
reslicing(FRETOut,resliced_directory,'Y',frequencyY,'FRET');
reslicing(FRETOut,resliced_directory,'Z',frequencyZ,'FRET');

 

end