The values of a frequency domain function represent how much of that frequency is "in" the function. For example, if you would take the fourier transform of a sine wave, you would get a delta function in the frequency domain: there's a lot of some specific frequency in that function.
Now, this is quite a simple way of saying it; we can deduct quite a bit more from the value, such as the "phase" of that frequency component. But it's also a bit more difficult, since for many signals the fourier transform is not composed of delta functions but a continuous function.
Don't worry too much about it, it takes a while getting used to the idea of fourier transforms. It's only now in my third year in Electrical Engineering that they really feel natural, I must say.
Here are short observations for your question. I think I'll add detailed explanation later.
Why frequency graphs are square?
This comes from nature of multidimensional DFT. One-dimensional DFT is just an expansion of vector in basis of complex harmonics. Similar construction can be applied for 2D images, but here we should take something like a tensor product of basis 1D signals along "$x$" and "$y$" axis of matrices. That'll be the basis for the space of 2D images. If we have a vector of length $N$, then we have $N$ vectors in Fourier basis, each of them stands for different frequency. When we have a $N \times N$ image and we take a tensor product of basic signals, we'll have $N^2$ combinations of frequencies along "$x$" and "$y$". Each vector from the Fourier basis for an image can be determined by its frequencies, which means form the tuple $(j, k)$. Frequency image shows the absolute value of expansion coefficient before the vector with frequencies $(j, k)$. So, that's why it's a square; generally it should be rectangular.
Why are the graphs point symmetrical around the center?
It's the consequence of formulas that used for Fourier expansion. For 1D DFT, coefficients before $f_j$ and $f_{N-j}$ are such that they have the same absolute value (this is true for dot products with real vector, as shown further). This property is inherited by 2D DFT and you see that $(j, k)$, $(N-j, k)$, $(N-j, N-k)$ and $(j, N-k)$ have the same absolute value.
In image 3, why is it a dense grid with 2 bright lines?
This image has structure that is much closer to a tensor product of 1D signals than the previous. And that 1D signal can be represented with very few harmonics. It would be interesting to see, how size of square affects the number of harmonics with big coefficients.
ADDITION
Here I'll try to explain rigorously some things that vere mentioned earlier.
I don't know how much you know about DFT so I'll try to cover as much as possible.
So, if we have $N$-dimensional comlex vector space, the Fourier basis can be given by following formula:
$$ F_m (n) = \frac{1}{\sqrt{N}} \exp{\frac{2\pi i \cdot m \cdot n}{N}},$$
where $m$ stands for the $m$-th basis vector and $n$ is a component index (both are integer and $0 \leqslant m,\; n \leqslant N-1$).
It is well known that Fourier basis is also an orthonormal basis with respect to standard dot product in $\mathbb{C}^n$, i.e. $\langle u, v \rangle = \sum \limits_{i} u_i \overline{v_i}$. So, then we are looking for expansion coefficients we can use dot product to find them. Here comes detailed explanation for second question. If we took a real vector $v$ and compare it's scalar product with $F_j$ and $F_{N-j}$:
$$ \langle v, f_j \rangle = \sum\limits_{p=0}^{N-1} v_p \overline{f_j(p)} =
\sum\limits_{p=0}^{N-1} v_p \overline{\frac{1}{\sqrt{N}} \exp{\frac{2\pi i \cdot j \cdot p}{N}}} =
\frac{1}{\sqrt{N}}\sum\limits_{p=0}^{N-1} v_p \exp{\left ( -\frac{2\pi i \cdot j \cdot p}{N} \right ) }$$
$$ \langle v, f_{N-j} \rangle = \sum\limits_{p=0}^{N-1} v_p \overline{f_{N-j}(p)} =
\sum\limits_{p=0}^{N-1} v_p \overline{\frac{1}{\sqrt{N}} \exp{\frac{2\pi i \cdot (N-j) \cdot p}{N}}} =
\frac{1}{\sqrt{N}}\sum\limits_{p=0}^{N-1} v_p \exp{\left ( \frac{2\pi i \cdot j \cdot p}{N} \right ) } $$
Since all coefficients of $v$ are real numbers we see that $\langle v, f_{N-j} \rangle = \overline{\langle v, f_{j} \rangle}$ and these dot products have the same absolute value.
Now I'll explain what I meant by tensor product of basis vectors in case of 2D images. By tensor product $u \otimes v$ I meant matrice $A$ which coefficient $A_{i, j}$ is equal to $u_i v_j$. Now we have two indices to locate the element of image and two indices for the frequencies. The formula for $(m, n)$ element of $(j, k)$ Fourier basis vector is following:
$$ F_{j,k} (m,n) = \frac{1}{\sqrt{N_1}} \exp{\frac{2\pi i \cdot j \cdot m}{N_1}}
\cdot \frac{1}{\sqrt{N_2}} \exp{\frac{2\pi i \cdot k \cdot n}{N_2}},$$
if we deal with $N_1 \times N_2$ matrices, or
$$ F_{j,k} (m,n) = \frac{1}{N} \exp{\frac{2\pi i \cdot j \cdot m}{N}}
\cdot \exp{\frac{2\pi i \cdot k \cdot n}{N}} = F_j(m) \cdot F_k(n)$$
in case of $N \times N$ matrice.
We also can apply a standard dot product to a space of 2D images if we put
$$ \langle u, v \rangle = \sum\limits_{m=0}^{N-1} \sum\limits_{n=0}^{N-1} u(m, n) \overline{v(m,n)}.$$ This nicely agrees with embedding $\mathbb{C}^{N\times N}$ in $\mathbb{C}^{N^2}$. It also can be shown that $F_{k, j}$ also forms an orthonormal basis. Coefficients for expansion can be found with analogous formulas via dot product. As I said earlier the mentioned property for real images is inherited by 2D DFT.
Best Answer
Think of it as a spinning wheel, where a positive frequency means that the wheel is spinning counter clockwise (positive in a coordinate system) while a negative frequency means that the wheel is spinning the other way.
the wheel can be described as $e^{i \theta}$ and if $\theta$ increases, we move counterclockwise but if $\theta$ decreases or if we write $e^{-i \theta}$ it goes clockwise instead