Average GrayValues within irregular shapes using a shapefile

I am very new to this kind of stuff. To describe what I am looking for is:

1. Find all Pixels in a certain irregular shape with the help of an existing shape img.
2. Use all that gray values to find the average or median
3. Color all Pixels from 1. with that average value
4. Do it for all shapes.

So here is my input img:

INPUT

And here is the shape file:

Shape

What you see there is a multicrystalline crystal as infrared transmission photograph. Maybe somebody has already experience which that kind of stuff or can help me with some usefull links/snippets. I would prefer python language.

I dont want to start at the beginning, driving into the wrong direction :D I hope someone can help.