Function lap1_axis 

fn lap1_axis(
    f: &[f64],
    mask: &[u8],
    l: usize,
    a: usize,
    n_axis: usize,
    coord: usize,
    n_total: usize,
) -> f64

Compute second derivative along one axis using mask-adaptive stencil.

Matches MATLAB’s lap1_mex.c logic:

• idx = 2*G[l+a] + G[l-a] selects the stencil type: 3 = central, 2 = forward, 1 = backward, 0 = zero
• At domain boundaries: i=0 → forward, i=N-1 → backward

Arguments

• f - field values
• mask - binary mask
• l - linear index of current voxel
• a - stride for this axis (1 for x, nx for y, nx*ny for z)
• n_axis - total extent for this axis (nx for x, nxny for y, nxny*nz for z)
• coord - axis coordinate as linear offset (i for x, jnx for y, knx*ny for z)
• n_total - total number of voxels (only used for z boundary detection)