Okay, I read through this again and might understand
what you want. However it looks a bit shaky. Here
is my understanding, to see if we are on the same
page.

Given masks Ma and Mb for corresponding 3d+time datasets
V and W, the goal is to find mask values Ma_1,...,Ma_N,
Mb_1,...,Mb_M, such that if vector A = sum(Ma_i * Vi) over
the N voxels in mask A (and N vectors in V), and B is a
similar sum for mask Mb and 3d+time dataset W, then r(A,B)
is maximized. The solution would have N+M betas over the
2 masks.

But this seems identical to finding the least-squares
solution to the equation: 0 = sum(Ai*Vi) + sum(Bi*Wi).
In this light, we can see how to solve it in AFNI.
But we can also see trouble with the method.


To solve this...

Use 3dmaskdump to dump the voxel time series across the
masks to a text file (or 2 files). Then run a simple
linear regression to fit (-polort 0) this model to an
all-zero vector. The resulting beta weights represent
the solution.

The betas would have to be mapped back into the volume.


Why this seems a little ugly...

As the mask sizes grow, the correlation values would
rise "quickly". And the correlations should hit 1.0
once the total number of voxels across the 2 masks
equals the number of time points minus 1.

This seems very highly dependent on the mask sizes,
and could yield very high correlations for time series
that do not actually look alike (in terms of the ROI
averages or principle components, say).

I would be very cautious with this method.

- rick