Here's a procedure you should try in a first attempt :
You have RGB images, so each images is a 3*RC matrix, where R is the number of rows and C the number of columns.
We will consider each pixel of an image as a realization of a random vector $\mathbf{x}$ of dimention 3. Hence, having an image, you have RC realization of $\mathbf{x}$. With the classification exemple you gave, each classe has 4 images, so we have 4*R*C realization of $\mathbf{x}$. We are abble to estimate some stastical moments, that we hope will be good descriptor, such as mean, covariance etc ... (mean can be estimated with empirical mean or median).
Let's just stick with the mean for now. Let us call $m_1,m_2,m_3,m_4$ the means of each class you gave.
Now we have a new image we want to classify, we estimate it's 3-dimensional mean $\mathbf{\mu}$, then using nearest neighbour, the class is given by :
$$
C=\arg_i\min\left(\lVert\mathbf{\mu}-m_i \rVert _2\right)
$$
This is very basic.
You may improve this scheme using additional features like variance (diagonal of the covariance matrix) or covariance matrix. You may estimate variance or covariance using empirical variance or empirical covariance estimator (realy just the classical formulae). However the euclidian distance won't be a good distance since you will comparing inhomogenous vector (the features will be composed of mean + covariance), you'll want to use a more advanced distance such as Bhattacharyya distance for the special case of gaussian distribution.