function [f, VIDEO] = backproject2(q, rVec, phiVec, xVec, yVec, intpol)

[xGrid, yGrid] = meshgrid(xVec, yVec);
     f = zeros( size( xGrid ) );
    Nr = length( rVec );
  Nphi = length( phiVec );
deltaR = rVec( 2 ) - rVec( 1 );

fig = figure;

colormap gray(256);
switch intpol
    case 'nearest'
        
        for phiIx = 1 : Nphi
            pause(.1);
            phi = phiVec( phiIx );
            proj = q( :, phiIx );
            r = xGrid * cos( phi ) + yGrid * sin( phi );
            rIx = round( r / deltaR + (Nr + 1) / 2 );
            f = f + proj( rIx );
            
            imagesc( yVec, xVec, f );
            axis xy;
            axis image;
            colorbar;
            title( 'Reconstructed image' );
            
            VIDEO( phiIx ) = getframe;
        end
        
    case 'linear'
        
        for phiIx = 1 : Nphi
            pause(.1);
            phi = phiVec( phiIx );
            proj = q( :, phiIx );
            r = xGrid * cos( phi ) + yGrid * sin( phi );
            rIx = r / deltaR + (Nr + 1)/2;% Same as before but don't round 

            % for each index in rIx
            for x = 1 : size( rIx, 1 )
            for y = 1 : size( rIx, 2 )
                % Get the index
                p_id_now = rIx( x, y );

                % init the result
                p_now = 0;

                % If the index isn't round
                if ( p_id_now - round( p_id_now ) ~= 0 )
                    % get the ones before and after
                    p_id_before = floor( p_id_now );
                    p_id_after  =  ceil( p_id_now );

                    % Get the proj value before and after
                    p_before = proj( p_id_before );
                    p_after  = proj( p_id_after  );

                    % Distance before-after to normalize weights, alwayse 1 !
                    %p_id_dist = dist( p_id_before, p_id_after );

                    % Compute weight before and after as the distance to the
                    % current index
                    w_before = abs( p_id_now   - p_id_before );
                    w_after  = abs( p_id_after - p_id_now    );

                    % Interpolate the projection value
                    p_now = p_before * w_before + p_after * w_after;
                else
                    % If the index is round, just read it
                    p_now = proj( p_id_now );
                end

                % add the value in the back projection
                f( x, y ) = f( x, y ) + p_now;
            end
            end
            
            imagesc( yVec, xVec, f );
            axis xy;
            axis image;
            colorbar;
            title( 'Reconstructed image' );
            
            VIDEO( phiIx ) = getframe;
        end
        
    otherwise
    error('Unknown interpolation.');
end

close( fig );

f = f * pi / (Nphi*deltaR);
