% SpatialIrsVideoMap_sps_image.m Spatial acoustic video mapping of an IR recording
% This script implements source acoustic mapping of IRs using array recordings 
% in sps format. Background is a static image. It saves video overlay with pseudocolor SPL maps.
%
% INPUT
% - sps_format.mat -> file containing the SPS format for IR, generated by c_generateEncodedFormat
% - panoramic_360_image.jpg
%
% OUTPUT
% acusticIrVideo.mp4
%
% Authors: Daniel Pinardi, Marco Binelli.
% University of Parma, 2017.
% Last review: 05/2023, Daniel Pinardi

clear variables
close all
clc

addpath( './lib' )                                                          % add to path the library with needed common files
a_configureFileNamesAndDirectories                                          % configure filenames, directories, resources

set( 0,'DefaultTextInterpreter', 'none' );                                  % setting for map titles
warning off                                                                 % to avoid warnings on image size


%% AUDIO PARAMETERS
% set fullscale value
% EM32_to_SPS32_2018    = 119
% EM32_to_SPS122_2018   = 122
fullScaleLevel          = 130;                                              % dB

% set EM32 gain to adjust FS level. Write -1 to avoid.
EM32gain                = -1;                                               % integer from 0 to 30

% Signal processing
trimEnable              = 1;                                                % enable trimming if 1
aWeightingEnable        = 1;                                                % enable to perform A filtering
timeSignalEnable        = 1;                                                % plot time signal under video maps. Disable for VR
nfft                    = 16384;                                            % reslution of the FFT
bufferSize              = 64;                                               % [samples]
stepSize                = 10;                                               % [samples]
autoselectOverlap       = 0;												% enable for VR video
frameRate               = 10;                                               % Default 30 (for VR) - min value for VLC is 6

% BAND PASS FILTER PARAMETERS -> 0 (one of the two) to disable
bandpassLowCutFreq      = 0;                                                % Low Cutoff Frequency [Hz]
bandpassHighCutFreq     = 2000;                                             % High Cutoff Frequency [Hz]


%% VIDEO/IMAGE PARAMETERS
% pointing correction: + -> to the right / - -> to the left
horizontalTilt          = 0; % degrees

% correction of the offset between optical and acoustic centers
imageOffsetCorrection   = 0;
radiusCorrection        = 1.5;      % [m]
verticalOffset          = 0.2;      % [m] 
outputWindowWidth       = 1024;                                             % desired width of background image

outputBlackWhite        = 1;                                                % 1 -> blackwhite, 0 -> color
alpha                   = 0.3;                                              % parameter for transparency mapping (0.4)

% Surf mapping
useFixedSplPlotBounds   = 0;                                                % if 0 each surface plot has color normalization 
fixedLowerLimit         = 80;                                               % user defined lower limit
fixedUpperLimit         = 90;                                               % user defined upper limit
noiseGate               = 20;                                               % max difference from maxlevel to plot
contourEnable           = 1;												% 1 to enable; 0 to disable. Show isolevel lines.
contourSpan             = 2; 												% [dB] Range of isolevel lines.   
valueTickEnable         = 0;												% 1 to enable; 0 to disable. Show values of isolevel lines.
tickSize                = 16;												% Font dimension.
tickBold                = 1;												% 1 to enable; 0 to disable. Font appearance.
colorbarEnable          = 0;												% 1 to enable; 0 to disable. Show colorbar
injectMetadata          = 0;												% 1 -> enable, 0 -> disable


%% LOAD DATA FROM FILES
% 32 or 122 SPS-irs
fprintf( '\nIrs elaboration started\n' );
filename                = encodedIrsMatFileName;                            % define filename to load and folder output
disp( [ 'Loading SPS format files: ' filename ] );
load( filename );
if exist( 'encodedFormat', 'var' ) == 1
    spsIrs              = encodedFormat;
    clear encodedFormat
end
virtualMicNumber        = size( spsIrs, 2 ); 
if virtualMicNumber == 32
    prompt              = 'SPS type (0 -> capsule directions ; 1 -> balloon): ';
    sps32type           = input( prompt );
else
    sps32type           = 2;
end

% Background 360° image 
disp( ['Loading background image file: ' backgroundImageFile ] );
backgroundImage         = imread( backgroundImageFile );


%% IMAGE PROCESSING
% fix width taking into account screen resolution scaling 
videoWindow				= figure( 'units', 'pixel', 'innerposition', [ 0 0 1000 500 ], 'menubar', 'none' );
F						= getframe( videoWindow );
widthsize				= size( F.cdata );
close( videoWindow );
if widthsize(2) ~= 1000
   ScreenScaleFactor	= widthsize( 2 ) / 1000;
   outputWindowWidth	= outputWindowWidth / ScreenScaleFactor;
end

% check monitor dimension and resize if necessary
[ outputWindowWidth ]   = screenSize( outputWindowWidth, timeSignalEnable );  

% check image aspect ratio and perform correction if needed ( output format is 2x1 ) and resize image
[ backgroundImage, backgroundHeight, backgroundWidth, resizeFactor ] = image2x1resize( backgroundImage, outputWindowWidth );

% correction of offset between acoustic center and image center
if imageOffsetCorrection == 1
    fprintf( 'Perform image correction\n\n' );
    [ backgroundImage ] = imageVertOffsetCorr( backgroundImage, radiusCorrection, verticalOffset );
end

% if outputBlackWhite == 1
%     background          = (1+alpha) * flipud( backgroundImage( :, :, [ 1 1 1 ] ) ); % flip and conversion to black white
% else
%     background          = (1+alpha) * flipud( backgroundImage );                    % flip
% end

if outputBlackWhite == 1
    background          = flipud( backgroundImage( :, :, [ 1 1 1 ] ) ); % flip and conversion to black white
else
    background          = flipud( backgroundImage );                    % flip
end


%% A-weighting
if aWeightingEnable == 1
    fprintf( 'Performing SPS format A-weighting.\n' );
    AWeighting          = weightingFilter( 'A-weighting', fs );             % use Matlab built-in A_weigth filter
    for i = 1 : virtualMicNumber
        spsIrs( :, i ) = AWeighting( spsIrs( :, i ) );
    end
end


%% TRIMMING
spsMean                 = mean( spsIrs, 2 );                                 % mean signal for plot

if trimEnable == 1
    
    fig                 = figure( 'windowstate', 'maximized' );
    hax                 = axes;
    plot( spsMean( :, 1 ) )
    xlim( [ 0 size( spsMean, 1 ) ] )
    grid on
    xlabel( 'Samples' )
    title( 'Audio recorded' )

    prompt              = 'AUDIO samples to skip ( 0 is accepted ): ';
    noiseStart          = input( prompt );
    line( [ noiseStart noiseStart ], get( hax, 'YLim' ), 'Color', [ 1 0 0 ] )
    prompt              = 'last AUDIO sample: ';
    noiseEnd            = input( prompt );

    close all

	spsIrs				= spsIrs( 1 + noiseStart : noiseEnd, : );      % trim noise
    
else % no trimming
    
    fprintf( 'Continuing without trimming.\n' );
	figure( 'windowstate', 'maximized' ) 
    plot( spsMean );
    xlim( [ 0 size( spsMean, 1 ) ] )
    xlabel( 'Samples' )
    grid on
    title( 'Untrimmed signal' );
    fprintf( 'Check signal...then press any key to continue or Ctrl-C to abort\n\n' );
    pause
    
end
close all

if autoselectOverlap == 1
    stepSize            = floor( fs / frameRate );							% overwrite
    overlap             = ( 1 - stepSize / bufferSize ) * 100;    
else
	overlap             = ( 1 - stepSize / bufferSize ) * 100;   
% 	stepSize            = floor( bufferSize * ( 1 - overlap ) );                % calculate stepSize
end

numbersFrame			= floor( length( spsIrs ) / stepSize );                 % integer number of frames
spsIrs					= spsIrs( 1 : stepSize * numbersFrame, : );             % update variable
fprintf( [ 'Overlap is ' num2str( overlap ) '%%.\n' ] )
fprintf( [ 'Frames are ' num2str( numbersFrame ) ' \n' ] )
fprintf( [ 'Video will be ' num2str( numbersFrame / frameRate ) 's \n' ] )


%% FULL SCALE LEVEL
if EM32gain >= 0 && EM32gain <= 30
	fullScaleLevel		= 0.0004*( EM32gain^3 ) - 0.0078*( EM32gain^2 ) - 0.8821*( EM32gain ) + fullScaleLevel; % dB
end

%% PADDING
zeroPadMatrix			= zeros( bufferSize/2, virtualMicNumber );              % center of buffer must be at 0
spsIrs					= [ zeroPadMatrix; spsIrs; zeroPadMatrix ];             % apply zero padding
spsMean					= mean( spsIrs, 2 );                                    % update value


%% Build grid
surfpoint				= 360;
if virtualMicNumber == 32 && sps32type == 0
    fprintf( 'Building 32 virtual mics grid - EM directions\n\n' );
    [ extX, extY, extXI, extYI, XI, YI, X, Y ] = GridBuilderEM32( backgroundWidth, backgroundHeight, surfpoint, virtualMicNumber );
elseif virtualMicNumber == 32 && sps32type == 1
    fprintf( 'Building 32 virtual mics grid\n\n' );
    [ extX, extY, extXI, extYI, XI, YI, X, Y ] = GridBuilder32( backgroundWidth, backgroundHeight, surfpoint, virtualMicNumber );
elseif virtualMicNumber == 64
	fprintf( 'Building 64 virtual mics grid\n\n' );
	[ extX, extY, extXI, extYI, XI, YI, X, Y ] = GridBuilder64( backgroundWidth, backgroundHeight, surfpoint, virtualMicNumber );
elseif virtualMicNumber == 122
    fprintf( 'Building 122 virtual mics grid\n\n' );
    [ extX, extY, extXI, extYI, XI, YI, X, Y ] = GridBuilder122( backgroundWidth, backgroundHeight, surfpoint, virtualMicNumber );
end



%% DATA PROCESSING
fprintf( 'Data processing started\n\n' );
if nfft < bufferSize
    nfft				= bufferSize;
end

freq					= linspace( 0, fs/2, nfft/2 + 1 );
tmp						= abs( freq - bandpassLowCutFreq );
[ ~, idxLowBP ]			= min( tmp );                                           % index of closest value
tmp						= abs( freq - bandpassHighCutFreq );
[ ~, idxHighBP ]		= min( tmp );
freqWindowing			= [ zeros( virtualMicNumber, idxLowBP ) ...
							ones( virtualMicNumber, idxHighBP - idxLowBP + 1 ) ...
							zeros( virtualMicNumber, ( nfft / 2 + 1 ) - idxHighBP - 1 ) ]';

if bandpassHighCutFreq == 0 || bandpassLowCutFreq == 0
    freqWindowing( freqWindowing == 0 ) = 1;
end
 
fftNormFactor			= sqrt( nfft * bufferSize );
splMatrix				= zeros( surfpoint/2, surfpoint, numbersFrame );        % preallocation

% Percentage initializing
oldCompletedPercentage  = -1;
fprintf( '\nFrames processing: ' )
varStr                  = '0%%';
varStrLen               = length( varStr );
fprintf( '%s', varStr )

for blockIndex = 1 : numbersFrame
  
    %%%%%%%%%%% PERCENTAGE PROCESSING %%%%%%%%%%%  
    completedPercentage = round( blockIndex/numbersFrame*100 );
    if oldCompletedPercentage ~= completedPercentage
        for i = 1 : varStrLen
            fprintf( '\b' )
        end
        varStr          = sprintf( '%d%%', completedPercentage );
        varStrLen       = length( varStr );
        fprintf( '%s', varStr )
    end
    oldCompletedPercentage = completedPercentage;
    
    % perform FFT on a chunk and compute RMS with complex sum
    frfChunk			= fft( spsIrs( ( blockIndex-1 ) * stepSize + 1 : ...
							  ( blockIndex-1 ) * stepSize + bufferSize, : ), nfft, 1 ) / fftNormFactor;
    frf					= 2 * frfChunk( 1 : nfft/2 + 1, : ) .* freqWindowing;

    complexSum			= real( frf ).^2 + imag( frf ).^2;
    rmsChunk			= ( sum( complexSum, 1 ) )' ./ 2;
    splChunk			= ( 10 * log10( rmsChunk ) + fullScaleLevel )';
    
    if virtualMicNumber == 122 || ( virtualMicNumber == 32 && sps32type == 1 )
        [ splChunk ]	= extendSPLdir( splChunk ); % replicate SPL value on poles
	elseif virtualMicNumber == 64
		[ splChunk ]	= extendSPLdir64( splChunk );
    end
    
    % interpolation
    [ ZI ] = splInterpolation( splChunk, virtualMicNumber, extX, extY, extXI, extYI, surfpoint, sps32type );
    
    splMatrix( :, :, blockIndex )   = ZI;

end

% pointing correction
correctionIndex			= round( horizontalTilt * surfpoint / 360 );
if correctionIndex > 0  
    splMatrix			= splMatrix( :, [ correctionIndex+1 : surfpoint, 1:correctionIndex ], : );
elseif correctionIndex < 0
    splMatrix			= splMatrix( :, [ surfpoint+correctionIndex+1 : surfpoint, 1:surfpoint+correctionIndex ], : );
end

upperLimit				= max( max( max( splMatrix ) ) );                       % find absolute max level
lowerLimit				= upperLimit - noiseGate;                               % range to plot

if useFixedSplPlotBounds == 1
    v					= fixedLowerLimit : contourSpan : fixedUpperLimit;      % iso-line levels
    splMatrix( splMatrix < fixedLowerLimit ) = 0;                           % set to 0 SPL values < lowerLimit
else                         
    v					= lowerLimit+1 : contourSpan : upperLimit;              % iso-line levels
    splMatrix( splMatrix < lowerLimit ) = 0;                                % set to 0 SPL values < lowerLimit
end
fprintf( 'Data processing ended correctly.\n\n' )


%% VIDEO PROCESSING
fprintf( 'Initializing video processing.\n\n' );
videoTitle				= recordedFileName( 1:end-4 );                          % name of file video
videoFile				= VideoWriter( videoTitle, 'MPEG-4' );                  % initializes video
videoFile.FrameRate		= frameRate;                                            % set frame rate
videoFile.Quality		= 100;                                                  % set quality
open( videoFile );

cmap					= jet( 64 );
cmap( 1, : )			= [ 0 0 0 ]; % black full scale

% Percentage initializing
oldCompletedPercentage  = -1;
fprintf( '\nFrames processing: ' )
varStr                  = '0%%';
varStrLen               = length( varStr );
fprintf( '%s', varStr )

if timeSignalEnable == 1
    
    videoHeigth			= outputWindowWidth / ( 2 * 0.8 );                      % define video window
    videoWindow			= figure( 'units', 'pixel', 'innerposition', [ 0 0 outputWindowWidth videoHeigth  ], 'menubar', 'none' );
    verticalLimit		= [ min( spsMean ) max( spsMean ) ];
    
    for blockIndex = 1 : numbersFrame
        
        %%%%%%%%%%% PERCENTAGE PROCESSING %%%%%%%%%%%  
        completedPercentage = round( blockIndex/numbersFrame*100 );
        if oldCompletedPercentage ~= completedPercentage
            for i = 1 : varStrLen
                fprintf( '\b' )
            end
            varStr		= sprintf( '%d%%', completedPercentage );
            varStrLen   = length( varStr );
            fprintf( '%s', varStr )
        end
        oldCompletedPercentage  = completedPercentage;
    
        ZI				= splMatrix( :, :, blockIndex );                        % extract current value to plot
        
        % subplot colormap
        plotBG			= subplot( 'Position', [ 0 0.2 1 0.8 ] );
        videoColorMap( background, backgroundWidth, backgroundHeight, XI, YI, ZI, lowerLimit, upperLimit, plotBG , ...
                       v, alpha, contourEnable, useFixedSplPlotBounds, fixedLowerLimit, fixedUpperLimit, valueTickEnable, ...
                       cmap, tickSize, tickBold, colorbarEnable );
        
        % subplot time signal
        plotTime		= subplot( 'Position', [ 0 0 1 0.2 ] );           
        timeSignalPlot( spsMean, verticalLimit, blockIndex, stepSize, bufferSize );

        writeVideo( videoFile, getframe( videoWindow ) );

    end
    
else
    
    videoHeigth			= outputWindowWidth / 2;                                % set video window dimension to 2:1 for VR
    videoWindow			= figure( 'units', 'pixel', 'innerposition', [ 0 0 outputWindowWidth videoHeigth  ], 'menubar', 'none' );

    for blockIndex = 1 : numbersFrame

        %%%%%%%%%%% PERCENTAGE PROCESSING %%%%%%%%%%%  
        completedPercentage = round( blockIndex/numbersFrame*100 );
        if oldCompletedPercentage ~= completedPercentage
            for i = 1 : varStrLen
                fprintf( '\b' )
            end
            varStr		= sprintf( '%d%%', completedPercentage );
            varStrLen   = length( varStr );
            fprintf( '%s', varStr )
        end
        oldCompletedPercentage = completedPercentage;
    
        ZI				= splMatrix( :, :, blockIndex );                        % extract current value to plot
    
        % subplot colormap
        plotBG			= subplot( 'Position', [ 0 0 1 1 ] );
        videoColorMap( background, backgroundWidth, backgroundHeight, XI, YI, ZI, lowerLimit, upperLimit, ...
                        plotBG, v, alpha, contourEnable, useFixedSplPlotBounds, fixedLowerLimit, fixedUpperLimit, ...
                        valueTickEnable, cmap, tickSize, tickBold, colorbarEnable );
        
        writeVideo( videoFile, getframe( videoWindow ) );

    end
end

close( videoFile );                                                         % finalizes video 
close all

if injectMetadata == 1
    source              = './lib/SpatialMediaMetadataInjector.exe';
    copyfile( source, './', 'f' )
    system( 'SpatialMediaMetadataInjector.exe' )
end
delete( 'SpatialMediaMetadataInjector.exe' )

rmpath( './lib' )
fprintf( '\n\nEnded correctly.\n\n' )

