EDIT: Note that this answer is completly wrong. See comments for proper answer.
You don't need sub2ind, simply use:
[rows, cols] = find(label == 41);
sub_feature_map = feature_map(rows, cols, :);
You can then reshape that into one column vector per pixel with:
sub_feature_map = reshape(permute(sub_feature_map, [3 2 1]), , size(sub_feature_map, 3));