function do_3dballgt( the_dir )
  
  if nargin < 1
    error( [ mfilename ' needs one input argument.' ] );
  end
    
  d = dir( fullfile( the_dir, 'seq*.ballgt' ) );
  
    % Parameters
  
  VIDEO_FRAMERATE = 25;

  % Interpolate between any two 3D measurements separated by less than this time interval
  MAX_INTERP_INTERVAL = 0.5; % seconds
  
  % To get some information about recordings: use the "seq_av163.m" script.
  CORPUS_PATH = '../..';
  addpath( CORPUS_PATH );

  % Loop through input files
  
  for i_file = 1:length( d )
    
    if ~d( i_file ).isdir  & ~isempty( strfind( d( i_file ).name, 'cam1' ) ) & isempty( strfind( d( i_file ).name, 'reprojected' ) ) & isempty( strfind( d( i_file ).name, 'interpolated' ) )
      
      IN_NAME1 = fullfile( the_dir, d( i_file ).name );
      IN_NAME2 = strrep( IN_NAME1, 'cam1', 'cam2' );
      IN_NAME3 = strrep( IN_NAME1, 'cam1', 'cam3' );
      
      if exist( IN_NAME2, 'file' ) & exist( IN_NAME3, 'file' )
	[ p,n,e ] = fileparts( d( i_file ).name );
	ind = strfind( n, 'cam1' );
	seqname = n( 1:ind(1)-2 );
	n = strrep( n, '-cam1', '' );

	% File that will contain 3D location estimates
	OUT_NAME3D = fullfile( the_dir, [ p n '.3dballgt' ] );
	
	% Files that will contain:
	% - original 2D manual measurements.
	% - between manual measurements: 2D reprojections of the 3D location estimates.
	%
	% The intent is to check visually whether or not the 3D reconstruction worked.
	OUT_NAME1 = strrep( IN_NAME1, '.ballgt', '-reprojected.ballgt' );
	OUT_NAME2 = strrep( IN_NAME2, '.ballgt', '-reprojected.ballgt' );
	OUT_NAME3 = strrep( IN_NAME3, '.ballgt', '-reprojected.ballgt' );
	
	disp( ' ' );
	disp( [ 'IN_NAME1:   "' IN_NAME1 '"' ] );
	disp( [ 'IN_NAME2:   "' IN_NAME2 '"' ] );
	disp( [ 'IN_NAME3:   "' IN_NAME3 '"' ] );
	disp( [ 'OUT_NAME3D: "' OUT_NAME3D '"' ] );
	disp( [ 'OUT_NAME1:  "' OUT_NAME1 '"' ] );
	disp( [ 'OUT_NAME2:  "' OUT_NAME2 '"' ] );
	disp( [ 'OUT_NAME3:  "' OUT_NAME3 '"' ] );
	
	% Load
	for a_cam = 1:3
	  ballgt = load( eval( sprintf( 'IN_NAME%d', a_cam ) ), '-ASCII' );
          % Compute 2D ball center on each camera
	  ballcentergt = [ ballgt( :, 1:3 ) mean( ballgt( :, [ 4 6 ] ), 2 ) mean( ballgt( :, [ 5 7 ] ), 2 ) ];
	  cam( a_cam ).ballgt = ballgt;
	  cam( a_cam ).ballcentergt = ballcentergt;
	end
	
	%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
	% ( 1 ) Load GT files and details (timecode)

	all_timestamps = [];
	
	o = seq_av163( seqname, CORPUS_PATH );
	
	if isempty( o )
	  disp( [ mfilename ': seq_av163() could not find sequence "' seqname '". Skipped.' ] );
	  break;
	end
	
	for a = 1:3
	  
	  % Save GT as is 

	  cam( a ).ballcentergt_orig = cam( a ).ballcentergt;
	  
	  % Convert frame # into time stamp
	  
	  cam( a ).ballcentergt( :, 1 ) = o.cam( a ).frame1_timestamp + ( cam( a ).ballcentergt( :, 1 ) - 1 ) / VIDEO_FRAMERATE;
	  
	  % Make sure the time stamp is unique
	  
	  cam( a ).ballcentergt( :, 1 ) = round( cam( a ).ballcentergt( :, 1 ) * 1e6 ) * 1e-6;
	  
	  % Keep only information about a visible mouth!
	  
	  cam( a ).ballcentergt = cam( a ).ballcentergt( find( cam( a ).ballcentergt( :, 3 ) == 1 ), : );
	  
	  % Fix the minor pixel variations in (X,Y) origin between sessions
	  
	  cam( a ).ballcentergt( :, 4 ) = cam( a ).ballcentergt( :, 4 ) - o.calib.cam( a ).delta_x_from_session08;
	  cam( a ).ballcentergt( :, 5 ) = cam( a ).ballcentergt( :, 5 ) - o.calib.cam( a ).delta_y_from_session08;

	  % Stack up the list of all timestamps
	  
	  all_timestamps = [ all_timestamps;   cam( a ).ballcentergt( :, 1 ) ];
	  
	end
	
	all_timestamps = sort( unique( all_timestamps ) );

	%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
	% ( 2 ) Reconstruct 3D mouth location when possible

	% Where we will store the result
	% - column 1: time stamp (in seconds)
	% - column 2: "is_manual" flag (0:interpolated or 1:manual)
	% - column 3 to 5: x,y,z coordinates
	
	ballcenter3dgt = zeros( 0, 5 );
	
	% Prepare things for "reconstruct_3d"
	for a = 1:3
	  cam( a ).Pmat    = o.calib.cam( a ).session08_Pmat;
	  cam( a ).K       = o.calib.cam( a ).session08_K;
	  cam( a ).kc      = o.calib.cam( a ).session08_kc;
	  cam( a ).alpha_c = o.calib.cam( a ).session08_alpha_c;
	end
	
	% Process data
	
	for t = all_timestamps(:).'
	  
	  % For each timestamp, load corresponding 2D measurements
	  
	  ind_cam = [ NaN NaN NaN ];
	  p2d = repmat( NaN, 3, 3 );
	  
	  is_manual = 1;
	  
	  for a = 1:3
	    ind = find( cam( a ).ballcentergt == t );
	    if length( ind ) == 1
	      ind_cam( a ) = ind;
	      p2d( :, a ) = [ cam( a ).ballcentergt( ind_cam( a ), 4:5 ).'; 1 ];
	      is_manual = is_manual & cam( a ).ballcentergt( ind_cam( a ), 2 );
	    end
	  end
	  
	  
	  % binary vector : 0 for unavailable measurement, 1 for available measurement
	  avail = ~isnan( ind_cam );
	  
	  which_avail = find( avail );
	  
	  p3d = [];
	  
	  if all( avail ) | all( avail == [ 1 1 0 ] ) | all( avail == [ 1 0 1 ] )
	    
	    % cams [ 1,2,3 ] or cams [ 1,2 ] or cams [ 1,3 ]
	    % NOT cams [ 2,3 ] because not precise enough in my opinion
	    
	    align_Pmat = eval( [ 'o.calib.rigid' sprintf( '%02d', which_avail ) '_Pmat' ] );
	    p3d = reconstruct_3d( reshape( p2d( :, which_avail ), [], 1 ), cam( which_avail ), align_Pmat );
	    
	    % Store the result
	    % the "is_manual" flag means that this 3D location estimate
	    % was derived from manual measurements
	    % on ALL cameras.
	    ballcenter3dgt( end+1, : ) = [ t is_manual p3d( 1:3 ).' ];
	    
	    
	  end
		  
	end 
	
	% Sort the matrix by increasing time stamp value
	[tmp, ind] = sort( ballcenter3dgt( :, 1 ) );
	ballcenter3dgt = ballcenter3dgt( ind, : );
	
	%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
	% ( 3 ) Reprojection on the image plane of each camera
	
	ind_notmanual = find( ~ballcenter3dgt( :,2 ) );
	m = ballcenter3dgt( ind_notmanual, 3:5 ).';
	m( 4,: ) = 1;
	p2d_notmanual    = project( m, cam, o.calib.rigid010203_Pmat );
	
	% Put the small pixel deviation (rel. to session08) back in
	for a = 1:3
	  p2d_notmanual( 1 + (a-1)*3, : ) = p2d_notmanual( 1 + (a-1)*3, : ) + o.calib.cam( a ).delta_x_from_session08;
	  p2d_notmanual( 2 + (a-1)*3, : ) = p2d_notmanual( 2 + (a-1)*3, : ) + o.calib.cam( a ).delta_y_from_session08;
	end

	t_new = ballcenter3dgt( ind_notmanual, 1 );
	
	% Insert the new points in the GT of each camera
	% Save it to a file
	
	for a = 1:3
	  
	  % Save GT as is
	  cam( a ).ballgt_orig = cam( a ).ballgt;
	  
	  % Translate manual gt back from "timestamp (in seconds)" to "frame number (integer)"
	  cam( a ).ballcentergt( :, 1 ) = 1 + round( ( cam( a ).ballcentergt( :, 1 ) - o.cam( a ).frame1_timestamp ) * VIDEO_FRAMERATE );
	  
	  % Add the new ( non manual ) measuremenets
	  frame_new = 1 + round( ( t_new - o.cam( a ).frame1_timestamp ) * VIDEO_FRAMERATE );
	  
	  % But only if they are not in the same frames as existing manual measurements 
	  % for THIS camera AND if they are not out of bound of this camera's recording
	  all_frames    = cam( a ).ballcentergt_orig( :, 1 );
	  ind_manual    = find( cam( a ).ballcentergt_orig( :, 2 ) );
	  frame_manual  = all_frames( ind_manual );
	  ind_ok = find( ~ismember( frame_new, frame_manual ) & ismember( frame_new, all_frames ) ); 
	  
	  n = length( ind_ok );
	  new_ballcentergt = [ frame_new(ind_ok(:)), repmat( [ 0 1 ], n, 1 ), p2d_notmanual( (1:2) + (a-1)*3, ind_ok ).' ];
	  
	  % Concatenate old values that we keep with new value that we add
	  
	  f1 = cam( a ).ballcentergt_orig( :,1 );
	  f2 = new_ballcentergt( :,1 );
	  
	  ind1 = find( ~ismember( f1, f2 ) );
	  ind2 = find( ~ismember( f2, f1( ind1 ) ) );
	  cam( a ).ballcentergt = [ cam( a ).ballcentergt_orig( ind1, : ); new_ballcentergt( ind2, : ) ];
	  
	  % Sort by increasing time index
	  
	  [ tmp, ind ] = sort( cam( a ).ballcentergt( :,1 ) );
	  cam( a ).ballcentergt = cam( a ).ballcentergt( ind, : );
	  
	  % Output the 2D ball gt
	  ballgt = cam( a ).ballgt_orig;
	  
	  frame_orig = ballgt(:,1);
	  frame_new  = cam( a ).ballcentergt(:, 1 );
	  
	  % Sanity check
	  if length( frame_orig ) ~= length( frame_new )
	    error( [ mfilename ': inconsistent length of frame index vectors.' ] );
	  end
	  if ~all( frame_orig == frame_new )
	    error( [ mfilename ': inconsistent frame index vectors.' ] );
	  end
	  
	  dx = cam( a ).ballcentergt( :, 4 ) - cam( a ).ballcentergt_orig( :, 4 );
	  ballgt( :, 4 ) = ballgt( :, 4 ) + dx;
	  ballgt( :, 6 ) = ballgt( :, 6 ) + dx;
	  
	  dy = cam( a ).ballcentergt( :, 5 ) - cam( a ).ballcentergt_orig( :, 5 );
	  ballgt( :, 5 ) = ballgt( :, 5 ) + dy;
	  ballgt( :, 7 ) = ballgt( :, 7 ) + dy;
	  
	  cam( a ).ballgt = ballgt;
	  
	  % Store 2D reprojection error
	  ind = find( dx & dy );
	  cam( a ).repr_error = sqrt( mean( [ dx(ind).^2   dy(ind).^2 ] ) );
	  
	  disp( sprintf( '  cam%d: 3D->2D reprojection error: [ dx:%.3g, dy:%.3g ] pixels', a, cam( a ).repr_error ) );
	  
	end
	
	%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
	% ( 4 ) Save results

	save_3dballgt( OUT_NAME3D, ballcenter3dgt );
	for a_cam = 1:3
	  save_ballgt( eval( sprintf( 'OUT_NAME%d', a_cam ) ), cam( a_cam ).ballgt );
	end
	
      end  % If all 3 files exist
      
    end  % if ~d( i_file ).isdir  & ~isempty( strfind( d( i_file ).name, 'cam1' ) )
    
  end  % Loop through files