Calculate intersection (overlapping area) of two polygons (of which one a 35 sided irregular polygon) without the use of polybool function

I would like to calculate the coverage of the Pointer's gamut by any color space or display gamut in both the CIE 1931 xy chromaticity diagram and the CIE 1976 u'v' chromaticity diagram in Matlab R2013b. The only downside is that I can't use the polybool function of Matlab, because my student license does not include licensing for the Mapping Toolbox.

For examples see the images linked to below. The first is for CIE 1931 xy, the second for CIE 1976 u'v'.

http://tweakers.net/ext/f/l7WIFLy64EqLh8hgB9LkF1Tw/full.png

http://tweakers.net/ext/f/W4TVl7SVeoOQwsCB2bJRc2bg/full.png

The blue polygon is the spectral locus and line of purples (line segment connecting both ends of the curved part) and makes up the edge of the CIE 1931 xy and CIE 1976 u'v' chromaticity diagrams respectively.

The green polygon is the edge of the Pointer's gamut (or in this case an approximation thereof)

The red triangle in this case is the color space of the Rec. 2020 standard for 8k UHD / Super Hi-Vision, but can be any 3 to 6 sided convex irregular polygon as long as the vertices are inside or on the blue line.

What I would like to do is calculate the intersection of the red and green polygon, so that when I divide the result of that by the area of the green polygon I know the coverage of the Pointer's gamut. This seems fairly easy using the polybool function in Matlab, but without using that I have no idea where to start. Any ideas?

Probably not very helpful, but the code I used for the images is:

      X = xlsread('pointerdata.xlsx','specloc','B2:B323');
      Y = xlsread('pointerdata.xlsx','specloc','C2:C323');
      P = xlsread('pointerdata.xlsx','calculations','DF9:DF44');
      Q = xlsread('pointerdata.xlsx','calculations','DG9:DG44');
      L = xlsread('pointerdata.xlsx','calculations','DF47:DF50');
      M = xlsread('pointerdata.xlsx','calculations','DG47:DG50');
      plot(X,Y,P,Q,L,M)
      axis equal
      axis ([0 0.8 0 0.9])
      set(gca,'XTick',0:.1:.8)

The Excel-file is also attached.

[Matt J]

What Toolboxes do you have?

Answer 1

[Image Analyst]

Do you have the Image Processing Toolbox? If so, use poly2mask to create binary images of the images, then AND them together

b1 = poly2mask(x1, y1, rows, columns);
b2 = poly2mask(x2, y2, rows, columns);
overlap = b1 & b2;

[Kid]

Yes, I do have the Image Processing Toolbox. I changed the vertices of the red triangle to that of the Adobe RGB color space, so that it's easier to recognize the overlapping area from the original Pointer's gamut (the Rec. 2020 covers nearly all of it, so hard to recognize).

The result of that is http://tweakers.net/ext/f/8RodQhJDR7JfSwfWrGeKbArG/full.png

I then added this code to the original code:

    S = 1000*P;
    T = 1000*Q;
    J = 1000*L;
    K = 1000*M;
    b1 = poly2mask(S,T,900,800);
    b2 = poly2mask(J,K,900,800);
    overlap = b1 & b2;
    overlap = flipud(overlap);
    imshow(overlap)

Which results in http://tweakers.net/ext/f/2deJkwDhHL1g7x40dtOTC7MA/full.png

If you compare that to the other image you can see that that indeed is the intersection. So now the only step that's missing is calculating the area of the white part, how do I do that?

[Image Analyst]

See a related answer I gave: http://www.mathworks.com/matlabcentral/answers/113621#comment_192039

Answer 2

[Kid]

Fully solved now:

      X = xlsread('pointerdata.xlsx','specloc','B2:B323');
      Y = xlsread('pointerdata.xlsx','specloc','C2:C323');
      P = xlsread('pointerdata.xlsx','calculations','DF9:DF44');
      Q = xlsread('pointerdata.xlsx','calculations','DG9:DG44');
      L = xlsread('pointerdata.xlsx','calculations','DF47:DF50');
      M = xlsread('pointerdata.xlsx','calculations','DG47:DG50');
      format long
      A = polyarea(X,Y);
      B = polyarea(P,Q);
      S = 10000*P;
      T = 10000*Q;
      J = 10000*L;
      K = 10000*M;
      b1 = poly2mask(S,T,9000,8000);
      b2 = poly2mask(J,K,9000,8000);
      overlap = b1 & b2;
      overlap = flipud(overlap);
      b = bwboundaries(overlap);
      coverage = b{1,1};
      coverageX = coverage(1:end,2);
      coverageY = coverage(1:end,1);
      coverageX = coverageX / 10000;
      coverageY = (9000 - coverageY) / 10000;
      plot(X,Y,coverageX,coverageY)
      axis equal
      axis ([0 0.8 0 0.9])
      set(gca,'XTick',0:.1:.8)
      C = polyarea(coverageX,coverageY);
      A
      B
      C
      B/A
      C/B

[aosheng xuanyuan]

good job!