Possible indexing error in patch_grid() function

[MartinK84]

We are currently trying to understand the code and algorithm by looking at the Cartesian demo example. The example works fine but while going through the code we noticed something odd that could possibly be a bug.

We tried to visualize the patches returned by the patch_grid function by breaking into grappaPrep.m at line 132:
patch_cr = getPatch(gMDL.sMask, pSize, gMDL.pDist, 'grid');
What is odd that for patches located near the or after the center of k-space (right part of k-space) there are several patches which include a single voxel with a location that is not matching other patch indices.

For example, just picking index 10000 from patch_cr will return:

K>> patch_cr{10000,2}
ans =
  Columns 1 through 9
       10000       15597       15598       15599       15600       15797       15798       15799       15800
  Columns 10 through 17
       15997       15998       15999       16000       16197       16198       16199       16200

In this example, the indices 2 to 17 perfectly describe a patch with neighboring indices. However, the first index of 10000 does not match the patch at all.

The reason for the problem can be found in getPatch.m in function [patch_c] = patch_grid(sMask, pSize, p):
At line 49:
inds = [(1:size(ninds,1)).', ninds]; % add patch center indices into inds
At this point, the center of each individual patch is added to the other indices of each patch. It's that adding of the patch center that causes the single outlier voxels. The indices added there for some voxels seem to not match the assigned neighboring voxels from the ninds variable, particularly at or after reaching the center area in k-space where the alternating sample vs. non-sampled pattern is broken.

[tianrluo]

@MartinK84
This might have been a confusion.

The integer indices you list are, in fact, not direct indices to the input sMask.
I virtually separated unsampled and sampled points on the sMask, and indexed them separately.
Specifically, the arrays in patch_c(ii, 2) are formed in a format that, its first member indexes which unsampled point among all unsampled ones the patch is referring to as center;
Similarly, the rest members index which sampled points among all sampled ones the patch is referring to as neighbors.
You can check that, in that 200 * 200 sized k-space example, the first elements, indexing the centers, has a maximum of ~17000 across all patches.
That is far from 200 * 200 = 40000, just because that number is only indexing the unsampled points.

This separate indexing was designed intentionally during developments, such that at the reconstruction stage, I can have a proper sized array that only contain the reconstructed signals, see kPS0 in Grappa/grappa.m, where 0 was to indicate the signals inside are not sampled but reconstructed.
During developments, this design frees me from worrying whether the program is possibly replacing any actually sampled signal with a reconstructed one, which behavior is not wanted in GRAPPA.

[MartinK84]

Thank you very much for the clarification. I was also guessing that it was intentional but couldn't easily find where that first value is omitted when using the patches for weight calculation and when applying the weights. Will check again tomorrow with my student and go through the code again to see if we can spot where it's separated.

Thank you very much for sharing your code and answering the questions!

[tianrluo]

@MartinK84 You are welcome.

In grappa.m, line 38 might be where you are looking for the separation.
Then at line 44, only the sampled signals, kPS1 (1 for sampled, 0 o.w.), are passed into doRecon;
and line 94 does the indices separation.

Now I started to think there are indeed something about data format that are shared/coupled implicitly across my functions.
This reference implementation is not in a perfect shape in terms of decoupling, although I've put quite some energy in that way...

[MartinK84]

Thank you, [ctrInds, ngbInds] = deal(inds_c{ii}(:,1), inds_c{ii}(:,2:end)); is what I was looking for which clearly has it going from 2:end. And I guess since the calibration is using the shifts_c variable which does not contain the additional index as the first element nothing similar is required there.

Your reference implementation is not impossible to follow, but also not that easy. We will see how hard it will get once we mess around with the non-Cartesian part and how it goes once we try to throw in our own data. In a week or two or so ;-)

The reason why we started to look so deeply into the Cartesian implementation was to have an easier start for understanding the algorithm (which we now understand much better) and because we wanted to see if we can figure out a reason why there are still some rediual folding artifacts left in the Cartesian demo (more than I'm used to from other Grappa implementations).