Hi, Clément-

For calculating a meaningful Jacobian, I think you could do 2 things:
A) just take the NL warp, and get the Jacobian of what is there at each voxel. Simples.
B) remove the rigid body components of AFF (no rotation or translation), and combine the remaining shear+scale parts of AFF and combine those with NL.

Does A or B make more sense to you? (There might not be one right answer to this.) My thought for A is to to use just the "local", not global, parts of warping for the Jacobian estimate. My thought for B is that the shear+scale could be sensical to include in shape deformation comparisons. For A, you just have the final WARP dset to use---you could remove the average shift of each subbrick, to ensure that is not present in the data.

To get B, how would this sound:

# get approximate rotation part of affine file AFF, and invert it
cat_matvec   \
    -ONELINE \
    AFF -IP  \
    > AFF_IP

# combine the inverse rot matrix with the original one,
# to remove the rotation part of the original
cat_matvec   \
    -ONELINE \
    AFF      \
    AFF_IP   \
    > AFF_noR

# remove the translation part of what we just calculated
set aaa     = `cat AFF_noR`
set aaa[4]  = 0
set aaa[8]  = 0
set aaa[12] = 0
echo ${aaa} > AFF_noR_noT

# Combine that non-rot and non-translational affine mat with the NL warp
3dNwarpCat \
    -overwrite \
    -prefix NL_noR_noT \
    -warp1 NL \
    -warp2 AFF_noR_noT

... and then NL_noR_noT is useful for B's Jacobian estimation. NB: you might still want to remove any average delta_x, delta_y and delta_z that may be in the NL or NL_noR_noT part.

How does that seem?

--pt