Thursday, September 3, 2015

Simple Stretching: Mapping a circular disc to a square

We can map the circular disc to a square by just linearly stretching each point radially so that the rim of the circle matches the rim of the square.




The equations for the mapping are provided below

where sgn(x) is the signum function.

(u,v) are circular coordinates in the unit disc.
(x,y) are coordinates in the square [-1,1] x [-1,1]

[C++ implementation source code below]
----------------------------------------------------------------------------------

// sample code accompanying the paper:
// "Analytical Methods for Squaring the Disc"
// http://arxiv.org/abs/1509.06344
#include <stdio.h> #include <math.h> #define epsilon 0.0000001 inline float sgn(float input) { float output = 1.0f; if (input < 0.0) { output = -1.0f; } return output; } // Simple Stretching map // mapping a circular disc to a square region // input: (u,v) coordinates in the circle // output: (x,y) coordinates in the square void stretchDiscToSquare(float u, float v, float& x, float& y) { if ( (fabs(u) < epsilon) || (fabs(v) < epsilon)) { x = u; y = v; return; } float u2 = u * u; float v2 = v * v; float r = sqrt(u2 + v2); // a trick based on Dave Cline's idea // link Peter Shirley's blog if (u2 >= v2) { float sgnu = sgn(u); x = sgnu * r; y = sgnu * r * v / u; } else { float sgnv = sgn(v); x = sgnv * r * u / v; y = sgnv * r; } } // Simple Stretching map // mapping a square region to a circular disc // input: (x,y) coordinates in the square // output: (u,v) coordinates in the circle void stretchSquareToDisc(float x, float y, float& u, float& v) { if ( (fabs(x) < epsilon) || (fabs(y) < epsilon)) { u = x; v = y; return; } float x2 = x*x; float y2 = y*y; float hypothenusSquared = x2 + y2; // code can use fast reciprocal sqrt floating point trick // https://en.wikipedia.org/wiki/Fast_inverse_square_root
    float reciprocalHypothenus =  1.0f/sqrt(hypothenusSquared);
    
    float multiplier = 1.0f;
    // a trick based on Dave Cline's idea
    if (x2 > y2) {
        multiplier = sgn(x) * x * reciprocalHypothenus;
    } else {
        multiplier = sgn(y) * y * reciprocalHypothenus;
    }

    u = x * multiplier;
    v = y * multiplier;
}


int main()
{
    float x,y;
    float u,v;

    stretchSquareToDisc(0.0f,0.6f,u,v);
    stretchDiscToSquare(u,v,x,y);
    
    printf("%f %f\n",u,v);
    printf("%f %f\n",x,y);

    stretchDiscToSquare(0.857f, 0.514496f,x,y);
    stretchSquareToDisc(x,y,u,v);
    
    printf("%f %f\n",x,y);
    printf("%f %f\n",u,v);

    return 0;
}


7 comments:

  1. I did test with 100 values of x and y in range [-1,1]. I run from square to disc and disc to square for getting a set of y values back. I also run every y value in x values from square to disc. I found only some y values with some x values can be retrieved from disc to square. The total of pairs running in this test is 100 x 100. The total of distance error is so small, 2.708423763042589e-14, but will cost me in some how. So, can I conclude that because a simple stretching method is not equiareal as mentioned in your paper, there will be a negligible value of distortion hidding. Pls. comment of what I conclude

    ReplyDelete
    Replies
    1. The simple stretching map is not equareal! Use the Shirley-Chiu concentric map instead

      see:
      http://psgraphics.blogspot.com/2011/01/improved-code-for-concentric-map.html

      Delete
  2. Confused about data test. My y values in the range of x values, say [min, max] of y value is [0.009, 0.126] and of x value is [-1,1]. It is OK to run square to disc and disc to square. If yes, my x value can also be in range [-0.6,0.6] for [0.009, 0.126]. Pls. correct.

    ReplyDelete
    Replies
    1. All values inside the square have to be between [-1,1]. None of the mappings will work otherwise

      Delete
  3. from
    1. (fabs(x) < epsilon) || (fabs(y) < epsilon)) is when shape of squre is slightly changed to disc at all
    2. (fabs(u) < epsilon) || (fabs(v) < epsilon)) is when shape of disc is slightly changed to squre at all

    My question why you set epsilon value to 0.0000001? Pls. guide me.

    ReplyDelete
    Replies
    1. The epsilon value was set arbitrarily to avoid numerical imprecision. You can change the value of epsilon as long as it is a reasonably small number.

      Delete
  4. From AnonymousSeptember 23, 2018:
    It is like a majic faster than my technique. Using (fabs(x) < epsilon) || (fabs(y) < epsilon)) and (fabs(u) < epsilon) || (fabs(v) < epsilon)) are so work in both simple stretching mapping and FG-Squircular Mapping. Thx a lot for a guide.

    ReplyDelete