Sorry, I misunderstood what you wanted to do. You can still use the maskave command with the -mask SELF option, but first prepare a masked version of the data with a 3dcalc command or a modified version of the mask.

# mask data
3dcalc -a dset.nii -b mask.nii -prefix maskdset.nii -expr 'a*notzero(a)*notzero(b)'
3dmaskave -mask SELF -dump maskdset.nii

# modify mask to exclude zero voxels
3dcalc -a dset.nii -b mask.nii -prefix maskmod.nii -expr 'b*notzero(a)*notzero(b)'
3dmaskave -mask maskmod.nii -mrange 1 1 -dump dset.nii

Alternatively, you can use 3dROIstats with the -nzmean option to compute the mean not including zero voxels for each sub-brick in the dataset or the 3dBrickStat command with the -non-zero option. Neither of these includes a dump of the individual voxel values though, but you can use a separate 3dmaskdump with combine mask and computed mask (-cmask) options.

# compute non-zero mean across ROI values and sub-bricks
3dROIstats -nzmean -mask mask.nii dset.nii

# compute single mean for all ROIs and sub-bricks
3dBrickStat -mask mask.nii -non-zero -mean dset.nii

# dump the voxel values in data inside mask that are non-zero in first sub-brick
3dmaskdump -mask mask.nii -cmask '-a dset.nii -expr a' dset.nii


I think these all will work, but you'll have to fiddle with the exact syntax for your data.