$\hat T_p$ is the complex fourier tranform. the bispectrum (triple correlation) is: \begin{equation} \hat T_{p,q}^{(3)}=\hat T_p.\hat T_q. \hat T_{-p-q} \end{equation}

If $\varphi$ is the phase of the complex fourier spectrum, $\beta$ is the phase of the complex bispectrum; The phase closure relation is:

$$ \exp(i \beta_{(p,q)})=\exp(i \varphi_p).\exp(i \varphi_q).\exp(i \varphi_{-p-q})$$ $$ \beta_{(p,q)}= \varphi_p + \varphi_q + \varphi_{-p-q} $$ $$ \varphi_{p+q}= \varphi_p + \varphi_q - \beta_{(p,q)} \ \ \ \ \ (1)$$

For real 1d signals, $\varphi_0=\varphi(0) = 0$, the imaginary part of the FT and the phase are odd hence:
$$\beta_{(0,1)}= \varphi_0 + \varphi_1 + \varphi_{0-1}$$ $$\beta_{(0,1)}= \varphi_1 - \varphi_{1}= 0 \ \ \ (?)$$

… maybe ok. But $\varphi_1 \ne 0$ and cannot be (or it's another signal…) Actually: $\beta_{(0,1)}=\beta_{(1,0)}=0$ is true. Same applies to: $\beta_{(n,0)} = \beta_{(0,n)}=0,\ \ \forall n$
Now, working on $\beta_{(n,1)}$:

$$\beta_{(1,1)}= \varphi_1 + \varphi_1 + \varphi_{-1-1}$$ $$\beta_{(1,1)}= 2 \varphi_1 - \varphi_2$$ $$\varphi_2 = 2 \varphi_1 - \beta_{(1,1)}$$
$$\beta_{(2,1)}= \varphi_2 + \varphi_1 + \varphi_{-2-1}$$ $$\beta_{(2,1)}= \varphi_2 + \varphi_1 - \varphi_{3}$$ $$\varphi_{3}= \varphi_2 + \varphi_1 - \beta_{(2,1)}$$ $$\varphi_{3}= 3 \varphi_1 - \beta_{(2,1)} - \beta_{(1,1)}$$ hence :
$$\varphi_{n}= n \varphi_1 - \beta_{(n-1,1)} - ... - \beta_{(1,1)}$$ $$\varphi_{n}= n \varphi_1 - \sum_{p=1}^{n-1}\beta_{(p,1)} \ \ \ \ \ n \gt 1 \ \ \ \ \ (2) $$ (1)

Same would apply working on $\beta_{(1,n)}$ instead of $\beta_{(n,1)}$, but that does not make independant determinations. Those will be coming from elsewhere.
Now, working on $\beta_{(n,2)}$:

$$\beta_{(3,2)}= \varphi_3 + \varphi_2 + \varphi_{-3-2}$$ $$\beta_{(3,2)}= \varphi_3 + \varphi_2 - \varphi_{5}$$ $$\varphi_{5}= \varphi_3 + \varphi_2 - \beta_{(3,2)} $$ avec: $$\varphi_{3}= \varphi_1 + \varphi_2 - \beta_{(1,2)}$$ $$\varphi_{5}= \varphi_1 + 2 \varphi_2 - \beta_{(3,2)}- \beta_{(1,2)} $$

Mais on n'ira pas tellement plus loin par là…

En pratique deux approches: si on chercher juste la phase sur une “ligne” la relation (2) suffit, mais on n'a aucune statistisque. C'est pratique en mise au point sans bruit, c'est une boucle. Sinon il faut appliquer la relation complète (1) en étant soigneux sur le support du calcul (les extrémités des boucles). Attempt to write (2) in 2d

For real 2d signals, the phase closure should be: $$ \beta(p_1,p_2,q_1,q_2)= \varphi(p_1,p_2) + \varphi(q_1,q_2) + \varphi(-p_1 -q_1 ,- p_2 -q_2) $$ $$ \varphi (p_1+q_1,p_2+q_2) = \varphi (p_1,p_2) + \varphi (q_1,q_2) - \beta(p_1,p_2,q_1,q_2) $$

For real 2d signals, $\varphi_{0,0}=\varphi(0,0) = 0$ (…), and the imaginary part of the FT and the phase are odd hence:

$$\beta(0,0,0,1)= \varphi(0,0) + \varphi(0,1) + \varphi(0,-1)$$ $$\beta(0,0,0,1)= \varphi(0,0) + \varphi(0,1) - \varphi(0,1) = 0$$ same for $$\beta(0,0,1,0)= \varphi(0,0) + \varphi(1,0) - \varphi(1,0) = 0$$ and also $\beta(0,1,0,0)=\beta(1,0,0,0)=0 $. (This is actually ok on my test image); From this, we cannot determine neither $\varphi(0,1)$ nor $\varphi(1,0)$. But if we know them already we can start the iteration at $$\varphi(1,1)=\varphi(0,1)+\varphi(1,0)-\beta(0,1,1,0)$$ (or symmetrical) beacuse the last term is not nul

$$\beta(0,1,0,1)= \varphi(0,1) + \varphi(0,1) + \varphi(0,-2)$$ $$\varphi(0,2)= 2 \varphi(0,1) - \beta(0,1,0,1)$$

$$\beta(0,1,0,2)= \varphi(0,1) + \varphi(0,2) + \varphi(0,-3)$$ $$\varphi(0,3)= \varphi(0,1) + \varphi(0,2) - \beta(0,1,0,2)$$ $$\varphi(0,3)= 3 \varphi(0,1) - \beta(0,1,0,1) - \beta(0,1,0,2)$$

hence the same kind of relationships:

$$\varphi(0,n)= n\ \varphi(0,1) - \sum_{p=1}^{n-1}\beta(0,1,0,p) = n\ \varphi(0,1) - \sum_{p=1}^{n-1}\beta(0,p,0,1) $$ $$\varphi(n,0)= n\ \varphi(1,0) - \sum_{p=1}^{n-1}\beta(1,0,p,0) = n\ \varphi(1,0) - \sum_{p=1}^{n-1}\beta(p,0,1,0)$$

But also:

$$\beta(0,1,x,1)= \varphi(0,1) + \varphi(x,1) + \varphi(-x,-2)$$ $$\varphi(x,2) = \varphi(0,1) + \varphi(x,1) - \beta(0,1,x,1)$$

$$\varphi(x,3) = \varphi(0,1) + \varphi(x,2) - \beta(0,1,x,2)$$ $$\varphi(x,3) = \varphi(0,1) + \varphi(0,1) + \varphi(x,1) - \beta(0,1,x,2)- \beta(0,1,x,1)$$ $$\varphi(x,3) = 2 \varphi(0,1)+ \varphi(x,1) - \beta(0,1,x,2)- \beta(0,1,x,1)$$ hence $$\varphi(x,n) = (n-1) \varphi(0,1) + \varphi(x,1) - \sum_{p=1}^{n-1}\beta(0,1,x,p)$$

The “quasi 1d” formula can be implemented right away and it works. Anyhow, it only allows to compute the (0,0) centered “cross” of phases. To compute the rest of it, the general formula must be used, and it is more a complication to stop the general algorithm to recompute that cross than to let it run freely. So in then end I removed this quasi 1d method and stay with the general formula. Unlike the 1d case, it is fundamental to initialize properly $\varphi(0,1)$ and $\varphi(1,0)$; otherwise nothing works.