How do I do phase unwrapping for a function in 2 variables for which I know the surface is increasing but I do not know the exact nature of the surface?

I have an unknown surface phi(i,k)on a 2 dimensional matrix such as 512x512 matrix. I know that it's minimum is phi(256,256)=0 and atan(tan(phi(256,256)))=0. it increases as radius from (256,256) increases. I know the values of atan(tan(phi(x,y), but they lie between -pi/2 to pi/2. How do I get the original continuously increasing surface phi(x,y)? I know I need to add pi whenever atan(tan(phi(x,y)) becomes less than 0 for the first time but not I am unable to figure out the code.