Alright, I've got a pretty good implementation in my branch https://github.com/WikiWatershed/mmw-geoprocessing/tree/tt/113/add-stac-cog-count-endpoint. It follows the recipe in Python / Jupyter here: https://gist.github.com/rajadain/1d3591a7a00c750466e3793bf1a9bcd2#file-lower-schuylkill-huc-10-test-2-ipynb.

The order of operations in Jupyter is:

Find tiffs in the STAC API that intersect the AoI
For each tiff, clip to AoI extent and reproject to EPSG:5070
- (really to the correct Albers projection based on where the AoI is, but for simplicity I'm using ConusAlbers EPSG:5070 here)
For each reprojected clip, mask to the AoI
Merge all masked clips
Calculate histogram using np.unique(data, return_counts=True)

This generates an output for the Lower Schuylkill HUC-10 like:

{
    "List(1)": 87250,
    "List(11)": 308877,
    "List(2)": 639920,
    "List(4)": 283,
    "List(5)": 173754,
    "List(7)": 2855276,
    "List(8)": 3126
}

The order of operations in GeoTrellis is:

Find tiffs in the STAC API that intersect with the AoI
Create a MosaicRasterSource from those tiffs, reprojected to EPSG:5070
Clip the RasterSource to the AoI's extent
Mask the Raster to the AoI
Calculate histogram using histogram.binCounts

This results in numbers that are quite close, but not exact:

{
    "List(1)": 105463,
    "List(11)": 390845,
    "List(2)": 730907,
    "List(4)": 525,
    "List(5)": 169751,
    "List(7)": 3591677,
    "List(8)": 5066
}