Hough transform

Description of the Hough transform

The input for the Hough transform are a set of points in a time frequency diagram that have been obtained from the demodulated FT.
The demodulation procedure has been performed with respect to a a certain sky location and certain spin-down parameters, in such a way that, if a source was located in the center of the patch and having the same spin-down parameters for which it has been demodulated, we will observe a set of points forming an horizontal line at the intrinsic frequency of the source f0. But because of the mismatch between the source parameters and the demodulated parameters we will observe a certain patterns in the time frequency diagram following the so called Hough transform master equation.

Suppose that the data has been demodulated for a given sky position N and spin-down parameters Fk, and imagine a source at another position inside the patch and with different spin-down parameters. The master equation tells us which will be the measured frequency \nu of the signal at a given time, that in general will be different from f0 , the intrinsic frequency of the signal.

The issue is to build histograms, i.e, a Hough map (HM), in the parameter space: for each intrinsic frequency f0, each residual spin-down parameter and each refined sky location inside the patch.
Notice, from the master equation, that the effect of the residual spin-down parameter is just a change in F0, and, at any given time, F₀ can be considered constant.
Also, the HM is an histogram, thus additive. It can be seen as the sum of several partial Hough maps constructed using just one periodogram.
Therefore, we can construct the HM for any f0 and spin-down value by adding together, at different times, PHM corresponding to different F₀ values.

In practice this means that in order to obtain the HM for a given frequency and all possible residual spin-down parameters, we have to construct a "CONE" of PHM around the frequency f₀. All of them coming from data demodulated with the same parameters.
The coordinates of the PHM locate the position of the source in the sky, and by summing along different directions owe refine the spin-down value.
When we want to analyze another frequency, for all possible spin-down parameters, we just need to add a new line to the cone (and remove another one) and then proceed making all the possible sums again.
For the case of only 1 spin-down parameter we have to sum following straight lines whose slope is related to the grid in the residual spin-down parameter. We can distinguish (at most) as many lines as the number of the different periodograms used.
We notice that the slope will be between -0.5 and +0.5. Therefore, independently of the line considered, we will be always adding together pairs of PHM at the same (horizontal) frequency.

So we can reduce the problem to this other picture that occupies less space in RAM.
Now the (new) lines have slopes between -1 and 1.
Since we would like to analyse different frequencies by shifting the "CONE", the expresssion of q(l,y) will become a bit more involved.
The code has to cope with the substitution of a line in the cone and this implies the use of pointers to the first element in each column as in a kind of circular buffer.

Given a time and f₀, implying a fixed value of F₀ and \xi, the master equation for a PHM can be expressed as: f=const1 + const2 * cos(phi), where: const1 = F0-\xi_x and const2 = |\xi| as it is used in the code. Thus, we can associate each value of cos(phi) to a frequency bin: fbin= "round"[(const1 + const2 * cos(phi))*Tcoh] where Tcoh is the inverse of the frequency resolution.
Going backt to the master equation, Let's see how to build one single PHM.
For a given value of F₀ and a measured frequency \nu (a selected peak), the source could be located anywhere on a circle with its center pointing in the same direction of \xi and characterized by the angle phi (the angle between n and \xi).
Since the Hough transform is performed on a set of spectra with discrete frequencies, if a peak on the spectrum appears at \nu it could correspond to any source with a demodulated frequency in a certain interval. As a consequence, the location of the sources compatible with F₀ and \nu is not a circle but an annulus with a certain width.
At this point we set the size of the pixel (sky location) to be smaller than half the width of the thinnest annulus.
An estimation of the average thickness of annuli tells us that the waste majority of annuli will be very thin, and that our algorithm should not be optimize for drawing thick annuli but thin ones. Also, the mapping to implement should be one with a uniform probability distribution in order to avoid discretization errors.
In order to remove border effects, we decided for a biunivoval mapping, meaning that a pixel can belong only to one annulus, (just touch by one peak of the spectrum). The criteria for the biunivocal mapping is that if and only if the center of the pixel is inside the annulus, then the pixel is enhanced.

Here there is an sketch of the zooming algorithm for the construction of a PHM.
Briefly, the zooming algorithm consists in: given a sky area (a patch or sub-patch), it checks if there are any annuli intersecting that area. There are three possible cases:

no anulli intersecting the patch (nothing has to be done)

all the patch is included in one annulus (enhance the number count in all pixels)

one or more annuli crossing the patch (the area needs to be refined: ZOOM)

For this pursopose, we need to calculate the minimum and maximum value of cos(phi) for that sky area. This can be very involved, and different subroutines have been developed in order to make the algorithm more efficient.
Depending upon the location of the center of the circle with respect to the patch, we distinguish(*):
corner case, inside case, opposite inside, parallel and meridian case.
(*further information to be provided upon request)

It is important that the number of pixels inside the patch to be a power of 2. We have been using a patch dimension of 256x256 pixels, that it was suitable for our choice of timescales.
Of course at this point we have already chosen a sky tiling to produce the PHM efficiently. It consists of changing coordinates so that the center of the
patch is located at (0,0) in alpha-delta (or in other coordinate system) and then we take parallels and meridians at constant space separation.

Last update by Alicia M. Sintes-Olives (26 Oct 2000)
http://www.aei-potsdam.mpg.de/~sintes/HOUGH/

	The input for the Hough transform are a set of points in a time frequency diagram that have been obtained from the demodulated FT. The demodulation procedure has been performed with respect to a a certain sky location and certain spin-down parameters, in such a way that, if a source was located in the center of the patch and having the same spin-down parameters for which it has been demodulated, we will observe a set of points forming an horizontal line at the intrinsic frequency of the source f0. But because of the mismatch between the source parameters and the demodulated parameters we will observe a certain patterns in the time frequency diagram following the so called Hough transform master equation.
	Suppose that the data has been demodulated for a given sky position N and spin-down parameters Fk, and imagine a source at another position inside the patch and with different spin-down parameters. The master equation tells us which will be the measured frequency \nu of the signal at a given time, that in general will be different from f0 , the intrinsic frequency of the signal.
	The issue is to build histograms, i.e, a Hough map (HM), in the parameter space: for each intrinsic frequency f0, each residual spin-down parameter and each refined sky location inside the patch. Notice, from the master equation, that the effect of the residual spin-down parameter is just a change in F0, and, at any given time, F₀ can be considered constant. Also, the HM is an histogram, thus additive. It can be seen as the sum of several partial Hough maps constructed using just one periodogram. Therefore, we can construct the HM for any f0 and spin-down value by adding together, at different times, PHM corresponding to different F₀ values.
	In practice this means that in order to obtain the HM for a given frequency and all possible residual spin-down parameters, we have to construct a "CONE" of PHM around the frequency f₀. All of them coming from data demodulated with the same parameters. The coordinates of the PHM locate the position of the source in the sky, and by summing along different directions owe refine the spin-down value. When we want to analyze another frequency, for all possible spin-down parameters, we just need to add a new line to the cone (and remove another one) and then proceed making all the possible sums again. For the case of only 1 spin-down parameter we have to sum following straight lines whose slope is related to the grid in the residual spin-down parameter. We can distinguish (at most) as many lines as the number of the different periodograms used. We notice that the slope will be between -0.5 and +0.5. Therefore, independently of the line considered, we will be always adding together pairs of PHM at the same (horizontal) frequency.
	So we can reduce the problem to this other picture that occupies less space in RAM. Now the (new) lines have slopes between -1 and 1. Since we would like to analyse different frequencies by shifting the "CONE", the expresssion of q(l,y) will become a bit more involved. The code has to cope with the substitution of a line in the cone and this implies the use of pointers to the first element in each column as in a kind of circular buffer.
Given a time and f₀, implying a fixed value of F₀ and \xi, the master equation for a PHM can be expressed as: *f=const1 + const2 cos(phi), where: const1 = F0-\xi_x and const2 = \|\xi\| as it is used in the code. Thus, we can associate each value of cos(phi)** to a frequency bin: fbin= "round"[(const1 + const2 * cos(phi))Tcoh] where Tcoh* is the inverse of the frequency resolution.	Going backt to the master equation, Let's see how to build one single PHM. For a given value of F₀ and a measured frequency \nu (a selected peak), the source could be located anywhere on a circle with its center pointing in the same direction of \xi and characterized by the angle phi (the angle between n and \xi). Since the Hough transform is performed on a set of spectra with discrete frequencies, if a peak on the spectrum appears at \nu it could correspond to any source with a demodulated frequency in a certain interval. As a consequence, the location of the sources compatible with F₀ and \nu is not a circle but an annulus with a certain width. At this point we set the size of the pixel (sky location) to be smaller than half the width of the thinnest annulus. An estimation of the average thickness of annuli tells us that the waste majority of annuli will be very thin, and that our algorithm should not be optimize for drawing thick annuli but thin ones. Also, the mapping to implement should be one with a uniform probability distribution in order to avoid discretization errors. In order to remove border effects, we decided for a biunivoval mapping, meaning that a pixel can belong only to one annulus, (just touch by one peak of the spectrum). The criteria for the biunivocal mapping is that if and only if the center of the pixel is inside the annulus, then the pixel is enhanced.
	Here there is an sketch of the zooming algorithm for the construction of a PHM. Briefly, the zooming algorithm consists in: given a sky area (a patch or sub-patch), it checks if there are any annuli intersecting that area. There are three possible cases: no anulli intersecting the patch (nothing has to be done) all the patch is included in one annulus (enhance the number count in all pixels) one or more annuli crossing the patch (the area needs to be refined: ZOOM) For this pursopose, we need to calculate the minimum and maximum value of cos(phi) for that sky area. This can be very involved, and different subroutines have been developed in order to make the algorithm more efficient. Depending upon the location of the center of the circle with respect to the patch, we distinguish(): corner case, inside case, opposite inside, parallel and meridian case. *(further information to be provided upon request)***
	It is important that the number of pixels inside the patch to be a power of 2. We have been using a patch dimension of 256x256 pixels, that it was suitable for our choice of timescales. Of course at this point we have already chosen a sky tiling to produce the PHM efficiently. It consists of changing coordinates so that the center of the patch is located at (0,0) in alpha-delta (or in other coordinate system) and then we take parallels and meridians at constant space separation.