A 4 × 4 grid of squares is filled in, with each of the 16 squares colored either black or white...

A 4 × 4 grid is filled in, with each of the 16 squares colored either black or white. Two colorings are regarded as identical if one can be converted to each other by performing any combination of flipping, rotating, or swapping the two colors (flipping all the black squares to white and vice versa). How many non-identical colorings are there?

I've figured out the number of invariances for each individual transformation but the combinations are a little confusing. Is there an easier way of solving this than just looking at each combination?

This is a case of Power Group Enumeration with the group permuting the slots being the eight symmetries $G_N$ of the $N\times N$ square and the group acting on the $Q$ colors being the symmetric group $S_Q$.

The cycle indices for $G_N$ were carefully documented and computed at the following MSE link I.

The cycle index of the symmetric group can be computed from the classical recurrence by Lovasz.

It then remains to apply the Power Group Enumeration formula / algorithm as documented at the following MSE link II.

We get for the case of coloring a square with at most two interchangeable colors $$1, 4, 51, 4324, 2105872, 4295327872, 35184441295872, \\ 1152921514807410688,\ldots$$ which is OEIS A182044 where a closed formula can be found.

For at most three colors we get the sequence $$1, 6, 490, 901012, 17653123147, 3126972103187548, \\ 4985402694248850150928, 71535079589434063488162675274, \\ 9238051838396620455005158025879826301,\ldots $$

Finally for at most four colors we get the sequence $$1, 7, 1555, 22396971, 5864091026091, 24595658827938966187, \\ 1650586719048786316922366635, 1772303994379887884341412962742479531, \\ 30447950777727144129269702965544605972525918891,\ldots$$

The sequence of colorings using any number of available interchangeable colors is also quite interesting and starts $$1, 7, 2966, 1310397193, \\ 579823814813639193, 477464341236883456112705749206, \\ 1340767144321669800049265230788088027597024910, \\ 21516767919669856366796245458194341840929552797762722429430679631, \ldots$$

An implementation of this algorithm is included below. The reader is invited to compute the closed formula from the algorithm specification which is not difficult to do but demands careful book-keeping.

with(combinat);

pet_cycleind_symm :=
proc(n)
option remember;

    if n=0 then return 1; fi;

    expand(1/n*add(a[l]*pet_cycleind_symm(n-l), l=1..n));
end;

pet_flatten_term :=
proc(varp)
local terml, d, cf, v;

    terml := [];

    cf := varp;
    for v in indets(varp) do
        d := degree(varp, v);
        terml := [op(terml), seq(v, k=1..d)];
        cf := cf/v^d;
    od;

    [cf, terml];
end;


pet_cycleind_grid :=
proc(N)
option remember;

    if type(N, even) then
        cind :=
        1/8*(a[1]^(N^2)+3*a[2]^(N^2/2)+
             2*a[1]^N*a[2]^((N^2-N)/2) + 2*a[4]^(N^2/4));
    else
        cind :=
        1/8*(a[1]^(N^2)+4*a[1]^N*a[2]^((N^2-N)/2)+
             a[1]*a[2]^((N^2-1)/2) + 2*a[1]*a[4]^((N^2-1)/4));
    fi;

    cind;
end;

gridcols :=
proc(N, Q)
option remember;
local idx_slots, idx_cols, res, a, b,
    flat_a, flat_b, cyc_a, cyc_b, len_a, len_b, p, q;

    if N > 1 then
        idx_slots := pet_cycleind_grid(N);
    else
        idx_slots := [a[1]];
    fi;

    if Q > 1 then
        idx_cols := pet_cycleind_symm(Q);
    else
        idx_cols := [a[1]];
    fi;


    res := 0;

    for a in idx_slots do
        flat_a := pet_flatten_term(a);

        for b in idx_cols do
            flat_b := pet_flatten_term(b);

            p := 1;
            for cyc_a in flat_a[2] do
                len_a := op(1, cyc_a);
                q := 0;

                for cyc_b in flat_b[2] do
                    len_b := op(1, cyc_b);

                    if len_a mod len_b = 0 then
                        q := q + len_b;
                    fi;
                od;

                p := p*q;
            od;

            res := res + p*flat_a[1]*flat_b[1];
        od;
    od;

    res;
end;

Addendum, Nov 2019. It is in fact not necessary to flatten the terms of the cycle indices and we may use them as is. This gives the following compact solution.

with(combinat);

pet_cycleind_symm :=
proc(n)
option remember;
local l;

    if n=0 then return 1; fi;

    expand(1/n*add(a[l]*pet_cycleind_symm(n-l), l=1..n));
end;

pet_cycleind_grid :=
proc(N)
option remember;
local cind;

    if type(N, even) then
        cind :=
        1/8*(a[1]^(N^2)+3*a[2]^(N^2/2)+
             2*a[1]^N*a[2]^((N^2-N)/2) + 2*a[4]^(N^2/4));
    else
        cind :=
        1/8*(a[1]^(N^2)+4*a[1]^N*a[2]^((N^2-N)/2)+
             a[1]*a[2]^((N^2-1)/2) + 2*a[1]*a[4]^((N^2-1)/4));
    fi;

    cind;
end;

gridcols :=
proc(N, Q)
option remember;
local idx_slots, idx_cols, res, term_a, term_b,
    v_a, v_b, inst_a, inst_b, len_a, len_b, p, q;


    if N > 1 then
        idx_slots := pet_cycleind_grid(N);
    else
        idx_slots := [a[1]];
    fi;

    if Q > 1 then
        idx_cols := pet_cycleind_symm(Q);
    else
        idx_cols := [a[1]];
    fi;

    res := 0;

    for term_a in idx_slots do
        for term_b in idx_cols do
            p := 1;

            for v_a in indets(term_a) do
                len_a := op(1, v_a);
                inst_a := degree(term_a, v_a);

                q := 0;

                for v_b in indets(term_b) do
                    len_b := op(1, v_b);
                    inst_b := degree(term_b, v_b);

                    if len_a mod len_b = 0 then
                        q := q + len_b*inst_b;
                    fi;
                od;

                p := p*q^inst_a;
            od;

            res := res +
            lcoeff(term_a)*lcoeff(term_b)*p;
        od;
    od;

    res;
end;