Use logical indexing to make one scatter() call each to plot all of the cases of each possibility. Then you can naturally legend() those.
Alternately, after you are done with the plotting, use plot(nan,nan) or line(nan,nan) once for each color and marker and fill style, and record the handles of each, and legend() that set of handles. Because of the nan coordinates those lines will not be drawn, but because they exist they can be legend()'d
The third possibility is that as you run through the plotting, each time that you detect that you have plotted a combination that you have not plotted before, record the handle of the scatter() -- so any one combination will be recorded only once. Then afterwards, legend() the recorded handles.
A fourth possibility is to record all of the scatter handles, and then afterwards post-process getting the attributes of each, and unique() on the combination of attributes with the two-output version of unique. The second output will give you indices into the handles to select to legend()
