The following Markov chain Monte Carlo algorithm would suffice, as it satisfies detailed balance and is ergodic. It would suffer from some sort of "critical slowing down" as the proportion of moves which are legal would decay as a power of the system size, but I don't have clear intuition about this.

I don't know if there's anything better in the literature. This is a "local" move; maybe there is a "global" move which is akin to the pivot move but which still gives a correct algorithm (e.g. rotating/reflecting/translating a piece of the polyomino and re-attaching it).

I'll just refer to the objects as "polyominoes", but for higher dimensions they could be referred to as $d$-dimensional polycubes.

1. Initialise system, perhaps starting with a "straight rod" of length $N$, starting at the origin and extending in the positive x-axis direction, or (better) starting with a d-dimensional cube with a few cells missing.
2. Choose a random occupied cell, and move it to a randomly chosen neighbouring position of the remaining $N-1$ cells. (For $Z^d$, there are $2d (N-1)$ options.) 
3. If the resulting configuration is a polyomino, accept the move, otherwise reject the move and keep the original configuration.
4. Return to step 2.

The algorithm generates translates of polyominoes, but that's trivial to partial out.

The Markov chain can generate any polyomino starting from a seed state, and is therefore ergodic. (Start with the rod seed state, and assume that the cell at the origin is the rightmost cell of an arbitrary polyomino that you want to construct. You can do this by taking the rightmost cell of the rod and building up the required configuration cell-by-cell to the left of the origin.)

The Markov chain satisfies detailed balance by inspection.

Therefore the Markov chain will sample polyominoes uniformly at random. It would take a bit of work to estimate the mixing time by running computer experiments and then coming up with a heuristic argument. I think it's unlikely that you will be able to get an exact result for the mixing time, but it may be possible to derive an upper bound.

One significant difficulty in implementation will be the efficient checking of connectivity for each attempted move. Naively, this could take CPU time $O(N)$, but it should be possible to do (much) better by using a dynamic bounding volume hierarchy (like my SAW-tree implementation of the pivot algorithm for self-avoiding walks) or a dynamic data structure which maintains information about connectivity (can't remember the name of the appropriate data structure for fully dynamic insertions and deletions of vertices, but they are related to link/cut trees).

Better algorithms are surely possible - e.g. the idea sketched above for a global move - but I think that this will suffice to get random polyominoes with say 1000 cells, even for naive connectivity checking. If you want to do this for a million cells then we might have to think a bit harder!