Everything @psadil said (and that I said in the linked post). One thing we did not address is motion censoring.
I’m not sure if this is the way everybody thinks about this, but sufficiently large motion produces artifacts that are induced by the resampling implicit in motion correction, and will not be fully eliminated by regressing out the motion parameters. In this case, you avoid these time points “infecting” your signal estimates by (in some way) removing them from your GLM.
In modeling hemodynamic responses to discrete events (as opposed to, e.g. rsfMRI), the strategy I’ve used is to zero-out the design matrix for the time points in question, and add a nuisance regressor (a column that is all zeros, except for a one at the time point in question) for each “censored” time point. In this case, I used ArtifactDetect to find time points with >(0.5 voxel-width) frame-wise motion or >3 sd frame-wise intensity jumps.
So to answer the question: Yes, I would use all three. Motion-correction is a baseline. For motion regression, note that the Friston 1996 paper recommends a 24-parameter autoregressive model (the DPARSF write-up contains a succinct description). And finally, you’ll need to choose a censoring strategy appropriate for your analysis, if there is one. The one I described above is not necessarily relevant to all analyses.