function interpolate_3dmouthgt( the_dir )

  if nargin < 1
    error( [ mfilename ' needs one input argument.' ] );
  end
  
  
  % Parameters
  VIDEO_FRAMERATE = 25;  % Hz
  MAX_INTERP_INTERVAL = 0.5;   % seconds
  XYZ_INTERP1_METHOD      = 'linear';  % does not matter much
  SPH_INTERP1_METHOD      = 'cubic';
  
  % To get some information about recordings: use the "seq_av163.m" script.
  CORPUS_PATH = '../..';
  addpath( CORPUS_PATH );
  
  file_list = dir( fullfile( the_dir, 'seq*.3dmouthgt' ) );
  
  for a = 1:length( file_list )
    
    if ~file_list(a).isdir & isempty( strfind( file_list(a).name, 'interpolated' ) ) & isempty( strfind( file_list( a ).name, 'reprojected' ) )
      
      in_mouth3d_filename = fullfile( the_dir, file_list(a).name );
     % Find the 3d ball gt file
      
      in_ball3d_filename = strrep( in_mouth3d_filename, 'mouth', 'ball' );
      if exist( in_ball3d_filename, 'file' )

	[ p,n,e ] = fileparts( in_mouth3d_filename );
	out_mouth3d_filename = fullfile( p, [ n '-interpolated' e ] );

	ind = strfind( n, '-person' );
	seqname = n( 1:ind-1 );
	personname = n( ind:end );
	out_mouth2d_name1 = fullfile( p, [ seqname '-cam1' personname '-interpolated-reprojected.mouthgt' ] );
	out_mouth2d_name2 = fullfile( p, [ seqname '-cam2' personname '-interpolated-reprojected.mouthgt' ] );
	out_mouth2d_name3 = fullfile( p, [ seqname '-cam3' personname '-interpolated-reprojected.mouthgt' ] );
	
	disp( ' ' );
	disp( [ 'in_mouth3d_filename:   "' in_mouth3d_filename '"' ] );
	disp( [ 'in_ball3d_filename:    "' in_ball3d_filename '"' ] );
	disp( [ 'out_mouth3d_filename:  "' out_mouth3d_filename '"' ] );
	disp( [ 'out_mouth2d_name1:     "' out_mouth2d_name1 '"' ] );
	disp( [ 'out_mouth2d_name2:     "' out_mouth2d_name2 '"' ] );
	disp( [ 'out_mouth2d_name3:     "' out_mouth2d_name3 '"' ] );
	
	mouthgt = load( in_mouth3d_filename, '-ASCII' );
	ballgt  = load( in_ball3d_filename, '-ASCII' );
	
	% Make sure the time is increasing
	
	[ tmp, ind ] = sort( mouthgt( :,1 ) );
	mouthgt = mouthgt( ind, : );

	[ tmp, ind ] = sort( ballgt( :,1 ) );
	ballgt = ballgt( ind, : );
	
	%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
	% ( 1 ) Mesure the distance between the mouth and the ball
	% This is just a quick check for the user to detect
	% if something is really wrong or not:
	% 
	% for most people this distance is about 20 to 25 centimeters
	% and the RMSE is very small.
	
	mouth.t = mouthgt( :, 1 );
	mouth.x = mouthgt( :, 3 );
	mouth.y = mouthgt( :, 4 );
	mouth.z = mouthgt( :, 5 );
	
	ball.t = ballgt( :, 1 );
	ball.x = ballgt( :, 3 );
	ball.y = ballgt( :, 4 );
	ball.z = ballgt( :, 5 );
	
	x = interp1( ball.t, ball.x, mouth.t, XYZ_INTERP1_METHOD );
	y = interp1( ball.t, ball.y, mouth.t, XYZ_INTERP1_METHOD );
	z = interp1( ball.t, ball.z, mouth.t, XYZ_INTERP1_METHOD );
	
	m_mouth = mouthgt( :, 3:5 );
	m_ball  = [ x(:) y(:) z(:) ];
	
	ind = find( ~isnan( x ) );

	% These three variables represent the knowledge we start with
	known.t = mouthgt( ind, 1 );
	known.mouth = m_mouth( ind, : );
	known.ball  = m_ball( ind, : );
	
	d = sqrt( sum( ( known.mouth - known.ball ) .^ 2, 2 ) );
	
	d_b_m = median( d );
	d_b_m_err = sqrt( mean( ( d - d_b_m ) .^ 2 ) );
	
	disp( sprintf( 'Distance ball-mouth:%.4g meters,   RMSE:%.4g', d_b_m, d_b_m_err ) );
	
	%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
	% ( 2 ) Interpolate in spherical coordinates,
	% taking the ball as the origin of the referent
	
	m = known.mouth - known.ball;
	[az,el,r] = cart2sph( m(:,1), m(:,2), m(:,3) );
	
	known.az = unwrap( az );
	known.el = unwrap( el );
	known.r  = r;

	unknown.t = [];
	ind = find( diff( known.t ) < MAX_INTERP_INTERVAL );
	for an_ind = ind(:).'
	  unknown.t = [ unknown.t  known.t( an_ind ):(1/VIDEO_FRAMERATE):known.t( an_ind+1 ) ];
	end
	unknown.t = sort( unique( [ unknown.t(:); known.t(:) ] ) );
		
	unknown.az = interp1( known.t, known.az, unknown.t, SPH_INTERP1_METHOD );
	unknown.el = interp1( known.t, known.el, unknown.t, SPH_INTERP1_METHOD );
	unknown.r  = interp1( known.t, known.r,  unknown.t, SPH_INTERP1_METHOD );
	
	%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
	% ( 3 ) Move back to the absolute referent
	% and cartesian coordinates
	
	[ dx, dy, dz ] = sph2cart( unknown.az, unknown.el, unknown.r );

	x = interp1( ballgt( :,1 ), ballgt( :, 3 ), unknown.t, XYZ_INTERP1_METHOD );
	y = interp1( ballgt( :,1 ), ballgt( :, 4 ), unknown.t, XYZ_INTERP1_METHOD );
	z = interp1( ballgt( :,1 ), ballgt( :, 5 ), unknown.t, XYZ_INTERP1_METHOD );
	
	unknown.ball  = [ x(:) y(:) z(:) ];
	unknown.mouth = unknown.ball + [ dx(:) dy(:) dz(:) ];
		
	out_mouthgt = [ unknown.t(:), 0 * unknown.t(:), unknown.mouth ];
	
	% Also put the old measures back in
	% ( we give priority to measurements directly derived from manual annotation )
	ind = find( ~ismember( out_mouthgt( :,1 ), mouthgt( :,1 ) ) );
	out_mouthgt = sortrows( [ out_mouthgt( ind, : );  mouthgt ] );
	
	% Save the result
	
	save_3dmouthgt( out_mouth3d_filename, out_mouthgt );
	
	%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
	% ( 4 ) Reproject on 2D image plane,
	% to check that we are not making gross errors
	
	cam = reproject_mouthgt3d( seqname, out_mouthgt, VIDEO_FRAMERATE, CORPUS_PATH );
	
	for a = 1:3
	  out_file = eval( sprintf( 'out_mouth2d_name%d', a ) );
	  base_file = strrep( out_file, '-interpolated-reprojected', '' );
	  
	  % Try to compare with manual 2D measurements,
	  % if we can still find them
	  
	  if exist( base_file, 'file' )
	  
	    % Load the manual measurements
	    base_mouthgt2d = sortrows( load( base_file, '-ASCII' ) );
	    ind_base = find( all( [ base_mouthgt2d(:, 2:3 )  ismember( base_mouthgt2d(:,1), cam( a ).mouthgt(:,1)) ], 2 ) );
	    f = base_mouthgt2d( ind_base, 1 );
	    ind_result = find( ismember( cam( a ).mouthgt(:,1), f ) );
	    
	    % Compute reprojection error
	    delta = base_mouthgt2d( ind_base, 4:5 ) - cam( a ).mouthgt( ind_result, 4:5 );
	    ind = find( all( delta, 2 ) );
	    disp( sprintf( 'cam%d: 3D->2D reprojection error: [ dx:%5.3f, dy:%5.3f ] pixels', a, sqrt( mean( delta( ind,:).^2 ) ) ) );

	    % Whenever possible put the manual measurements back in
	    % -> this way the "interpolated-reprojected" file can be updated
	    %    where necessary and re-fed to "s040827_3dmouthgt" and then to the present script
	    cam( a ).mouthgt( ind_result, : ) = base_mouthgt2d( ind_base, : );
	    
	  else
	    disp( sprintf( 'cam%d: 3D->2D reprojection error: [ could not find manual measurements ]', a ) );
	  end

	  % Save the result
		  
	  save_mouthgt( out_file, cam( a ).mouthgt );
	  

	  
	end
	
      end
      
    end
    
  end
  
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function cam = reproject_mouthgt3d( seqname, mouthgt3d, videoframerate, corpus_path )
  
  if nargin < 4
    error( [ mfilename ': reproject_mouthgt3d() needs four input arguments.' ] );
  end
  
  if size( mouthgt3d, 2 ) ~= 5
    error( [ mfilename ': reproject_mouthgt3d(): "mouthgt3d" must be a 5-column matrix.' ] );
  end
  
  % Load projection information
  o = seq_av163( seqname, corpus_path );
  
  p3d = mouthgt3d( :, 3:5 ).';
  p3d( 4,: ) = 1;
  
  for a = 1:3
    % Prepare things for "project"
    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  
  
  % Projection 3D -> 2D
  p2d = project( p3d, cam, o.calib.rigid010203_Pmat );
  
  % Put the small pixel deviation (rel. to session08) back in
  for a = 1:3
    p2d( 1 + (a-1)*3, : ) = p2d( 1 + (a-1)*3, : ) + o.calib.cam( a ).delta_x_from_session08;
    p2d( 2 + (a-1)*3, : ) = p2d( 2 + (a-1)*3, : ) + o.calib.cam( a ).delta_y_from_session08;
  end
  
  % Store the result
  t = mouthgt3d( :,1 );
  for a = 1:3
    % Frame index
    frame = 1 + round( ( t - o.cam( a ).frame1_timestamp ) * videoframerate );

    % "0" means is_manual  flag is set to 0
    % "1" means is_visible flag is set to 1
    z = 0 * frame;
    cam( a ).mouthgt = [ frame  z  1+z  p2d( (1:2) + (a-1)*3, : ).' ];
  end