You can use 3dQwarp with options "-maxlev 0 -noQ -allineate" to make it compute a warp with the affine part ("-allineate") plus a single 3D cubic polynomial ("-maxlev 0 -noQ"). These results should be similar to the mostly useless bilinear warp and will be available for use with other programs, such as 3dNwarpApply.
"-noQ" is a hidden option, which just means to skip the quintic polynomial steps at level 0 (the global patch). This option will make level 0 run pretty fast, since it will just do one cubic optimization.