The greatest common divisor is a well-defined notion in any unique factorization domain.
http://en.wikipedia.org/wiki/Greatest_common_divisor#The_gcd_in_commutative_rings
The greatest common divisor is unique up to multiplication by a unit; in the case of gaussian integers, this means one of ##\pm 1## or ##\pm i##.
A euclidean domain is a special case of a unique factorization domain, one in which the euclidean algorithm works, meaning essentially that we can divide any element ##a## by any nonzero element ##b## and get a quotient and a remainder, where the remainder is "smaller" than ##b##. Part of the definition of euclidean domain is a "norm" function which quantifies what we mean by smaller. It can be any function ##f## from the ring to the nonnegative integers, such that for every ##a## and every nonzero ##b##, there exist ##q,r## satisfying ##a = qb + r##, with either ##r = 0## or ##f(r) < f(b)##.
The gaussian integers are an example of a euclidean domain, with norm function ##f(x+iy) = x^2 + y^2##.
To address the OP's question, if ##a = 4 + 7i## and ##b = 1 + 3i## then for the first step of the euclidean algorithm, we seek ##q## and ##r## satisfying
$$4 + 7i = q(1+3i) + r$$
such that either ##r = 0## or the norm of ##r## is strictly smaller than ##1^2 + 3^2 = 10##. To find a candidate for ##q##, we can carry out the division as complex numbers (not constrained to be gaussian integers):
$$\frac{4+7i}{1+3i} = \frac{(4+7i)(1-3i)}{(1+3i)(1-3i)} = \frac{25 - 5i}{10} = \frac{5}{2} - \frac{1}{2}i$$
If we round the real part ##5/2## down to ##2## and the imaginary part ##-1/2## up to ##0##, we get a candidate ##q = 2 + 0i = 2##. Multiplying this candidate by ##1+3i##, we get ##2 + 6i##, and subtracting this from ##4+7i##, we end up with ##r = 2 + i##, which indeed satisfies the desired norm inequality ##2^2 + 1^2 < 1^2 + 3^2##.
Note that we could equally well have tried rounding the real part ##5/2## down to ##2## and the imaginary part ##-1/2## down to ##-1##, in which case we would have a candidate ##q = 2 - i##. Multiplying by ##1+3i## gives ##5 + 5i##, and subtracting this from ##4+7i##, we end up with ##r = -1 + 2i##, which also satisfies the norm inequality: ##1^2 + 2^2 < 1^2 + 3^2##. So this shows that the ##q## and ##r## at each step are not necessarily unique. This is OK, because at some point we must still end up with ##r=0## since at each step the norm of ##r## is a nonnegative integer which is strictly smaller than in the previous step.
Let's say we chose ##q=2## and ##r = 2+i## as in the OP. Then at the second step, we need to find ##q## and ##r## satisfying
$$1+3i = q(2+i) + r$$
where either ##r=0## or the norm of ##r## is strictly less than ##2^2 + 1^2 = 5##. To find a candidate ##q##, we can once again carry out the division in ##\mathbb{C}##:
$$\frac{1 + 3i}{2 + i} = 1 + i$$
No rounding is needed this time: we got an exact gaussian integer with no remainder. So this is the final step: ##q = 1 + i## and ##r = 0##, and the gcd is ##1+i## times any unit (##\pm 1## or ##\pm i##).