matlab 查找由一组直线创建的平面的所有分区

41zrol4v  于 2022-12-13  发布在  Matlab
关注(0)|答案(5)|浏览(180)

我有一组直线(形式为y = mx + B的线性函数)(共120条!),如果我把它们都画出来,它们就分割了R^2平面,这些直线不一定经过原点。
要找到由一组这样的行创建的所有分区,最有效的方法是什么?就我个人而言,我很难想出任何方法,更不用说有效的方法了。为了更清楚起见,我包括了以下只有4行的图像:x1c 0d1x
分区的一个例子是集合{(x,y)| -30x+28<= y && 60x+2 <= y <= 90x+7},它是由第一象限中的红、黄和绿色线创建的分区。另一个例子是{(x,y)|y <= -30x+28 && 5x+3 <= y <= 60x+2},它是第一象限中由蓝、红和绿线包围的三角形。
{(x,y)|5x+3 <= y <= -30x+28}是一个 non 分区的例子,它是由上面的绿色线和下面的蓝线围成的集合。这不是一个分区,因为它包含了几个分区(例如上面的第二个集合),或者与它重叠。然而,集合{(x,y)|5x+3 <= y <= -30x+28 && 90x+7 <= y}是一个分区。
所需的输出将是此类集合的完整列表:{(x,y)|y <= -30x+28 && 5x+3 <= y <= 60x+2},{(x,y)| -30x+28<= y && 60x+2 <= y <= 90x+7}...等。当然,它们不必以这种表示法给出。
我不知道如何解决这个问题,因此,不幸的是,我不能提供太多我所尝试的方法。理想情况下,我想用R,Python,Mathematica或MATLAB来做这件事,但在这一点上,我对任何选择都持开放态度。

**EDIT:**由于似乎有一个符号的问题,我将澄清一点。它足以简单地得到一个点的条件列表,这样所有满足该条件的点将精确地定义一个分区。例如,一个长的交集列表,将是很好的:y <= 5x+3 && y >= 90x+7 && y<= -30x+28是一个很好的定义分区的输出。当然,期望的输出是这样的分区的完整列表(如上所定义)。

hmae6n7t

hmae6n7t1#

求分区的数量遵循以下公式(当没有3条或更多的线相交于同一点时--这一假设贯穿于本文):

num_partitions = (n_lines * (n_lines + 1) / 2 ) + 1

一个解释是here,另一个解释是there
所需的输出将是此类集合的完整列表:{(x,y)|y <= -30x+28 && 5x+3 <= y <= 60x+2}, {(x,y)| -30x+28<= y && 60x+2 <= y <= 90x+7}...等...当然,它们不必用这种表示法给出。

这里缺少精确的符号是一个障碍。

下面是我的信封尝试的背面。如你所见,可以根据它们与每一行的相对位置来识别编号区域。
有5个空集,如果行顺序发生变化,则空集将不相同。

它可能会更容易划分平面上的一组点与线,试图确定哪个点属于哪个集;在这种情况下,探索2^n可能的分区并返回它们的内容将是很容易的。(比试图找到一个好的符号来标识抽象集要容易)

这并没有完全回答你的问题,但可能是一个很好的起点,有人能够/愿意进一步推动这一点。
这里有一个关于partitioning a set of points with two lines in the plane.的说明,这是一个不同的问题,但它的一些方法可能是有用的。

其他方法:

识别由线段形成的多边形,计算船体,确定点是否在该凸包中。

kd3sttzy

kd3sttzy2#

下面是Mathematica中的一个解决方案。该方法包括找到直线的交点、线段和分割,同时跟踪这些点连接到哪些直线上。

y[1] = 5 x + 3;
y[2] = 90 x + 7;
y[3] = -30 x + 28;
y[4] = 60 x + 2;

findpoints[n_] := Module[{},
  xp = DeleteCases[{{First[#1], First[#2]},
       Solve[Last[#1] == Last[#2], x]} & @@@
     DeleteCases[
      Tuples[Array[{#, y[#]} &, n], 2],
      {{x_, _}, {x_, _}}], {_, {}}];
  yp = y[#[[1, 1]]] /. #[[2, 1]] & /@ xp;
  MapThread[{#1, {#2, #3}} &,
   {xp[[All, 1]], xp[[All, 2, 1, 1, 2]], yp}]]

xyp = findpoints[4];

{xmin, xmax} = Through[{Min, Max}@
    xyp[[All, 2, 1]]] + {-0.7, 0.7};

outers = Flatten[Array[Function[n,
     MapThread[List[{{n, 0}, {##}}] &,
      {{xmin, xmax}, y[n] /.
    List /@ Thread[x -> {xmin, xmax}]}]], 4], 2];

xyp = Join[outers, xyp];

findlines[p_] := Module[{},
  pt = DeleteCases[
    Cases[xyp, {{p[[1, 1]], _}, _}], p];
  n = First@Cases[pt,
     {_, First@Nearest[Last /@ pt, Last[p]]}];
  {{First[p], First[n]}, {Last[p], Last[n]}}]

lines = Map[findlines, xyp];

(* boundary lines *)
{ymin, ymax} = Through[{Min, Max}@outers[[All, 2, 2]]];
{lbtm, rbtm, ltop, rtop} = {{xmin, ymin},
   {xmax, ymin}, {xmin, ymax}, {xmax, ymax}};
xminlines = Partition[Union@Join[{ymin, ymax},
      Cases[xyp, {_, {xmin, _}}][[All, 2, 2]]], 2, 1] /.
   x_Real :> {xmin, x};
xmaxlines = Partition[Union@Join[{ymin, ymax},
      Cases[xyp, {_, {xmax, _}}][[All, 2, 2]]], 2, 1] /. 
   x_Real :> {xmax, x};
lines2 = Join[Last /@ lines, xminlines, xmaxlines,
   {{lbtm, rbtm}}, {{ltop, rtop}}];

ListLinePlot[lines2]

(* add vertex points *)
xyp2 = Join[xyp, {
   {{SortBy[Cases[outers, {_, {xmin, _}}],
       Last][[-1, 1, 1]], -1}, ltop},
   {{SortBy[Cases[outers, {_, {xmax, _}}],
       Last][[-1, 1, 1]], -1}, rtop},
   {{SortBy[Cases[outers, {_, {xmin, _}}],
       Last][[1, 1, 1]], -1}, lbtm},
   {{SortBy[Cases[outers, {_, {xmax, _}}],
       Last][[1, 1, 1]], -1}, rbtm}}];

anglecalc[u_, v_] := Mod[(ArcTan @@ u) - (ArcTan @@ v), 2 π]

getlineangles[] := Module[{},
  (* find the angles from current line
     to all the linked lines *)
  angle = Map[
    anglecalc[{c, d} - {g, h}, # - {g, h}] &,
    union = DeleteCases[Union@Join[
        Last /@ Cases[lines2, {{g, h}, _}],
        First /@ Cases[lines2, {_, {g, h}}]],
      {c, d}]];
  Sort[Transpose[{N@angle, union}]]]

getpolygon[pt_, dir_] := Module[{},
  Clear[p];
  p[n = 1] = {{a, b}, {c, d}} = pt;
  (* find the angles from vector (0, -1) or (0, 1)
     to all the linked lines *)

  angle = Map[anglecalc[If[dir == 1, {0, -1}, {0, 1}], # - {c, d}] &,
    union = Union@Join[
       Last /@ Cases[lines2, {{c, d}, _}],
       First /@ Cases[lines2, {_, {c, d}}]]];
  lineangles = Sort[Transpose[{N@angle, union}]];
  (* next point *)
  p[++n] = {{e, f}, {g, h}} = First@
     Cases[xyp2, {_, lineangles[[1, 2]]}];

  While[Last[p[n]] != Last[p[1]],
   lineangles = getlineangles[];
   (* reset first point *)
   {{a, b}, {c, d}} = {{e, f}, {g, h}};
   (* next point *)
   p[++n] = {{e, f}, {g, h}} = First@
      Cases[xyp2, {_, lineangles[[1, 2]]}]];

  Array[p, n]]

len = Length[xyp];

polygons = Join[Array[(poly[#] = getpolygon[xyp[[#]], 1]) &, len],
   Array[(poly[# + len] = getpolygon[xyp[[#]], 2]) &, len]];

graphics = DeleteDuplicates /@ Array[Last /@ poly[#] &, 2 len];

sortedgraphics = Sort /@ graphics;

positions = Map[Position[sortedgraphics, #] &,
    DeleteDuplicates[sortedgraphics]][[All, 1, 1]];

unique = poly /@ positions;

poly2 = unique[[All, All, 2]];

poly2 = Delete[poly2,
   Array[If[Length[Intersection[poly2[[#]],
         Last /@ Take[xyp2, -4]]] == 4, {#}, Nothing] &,
    Length[poly2]]];

len2 = Length[poly2];

poly3 = Polygon /@ Rest /@ poly2;

Array[(centroid[#] = RegionCentroid[poly3[[#]]]) &, len2];

Show[Graphics[Array[{ColorData[24][#],
     poly3[[#]]} &, len2], AspectRatio -> 1/GoldenRatio],
 Graphics[Array[Text[#, centroid[#]] &, len2]]]

unique2 = Extract[unique,
   Flatten[Position[unique[[All, All, 2]], #] & /@ poly2, 1]];

makerelations[oneconnection_, areanumber_] := Module[{i},
  If[Intersection @@ oneconnection == {} ||
    (i = First[Intersection @@ oneconnection]) < 1,
   Nothing,
   centroidx = First[centroid[areanumber]];
   linepos = y[i] /. x -> centroidx;
   relation = If[linepos < Last[centroid[areanumber]],
     " >= ", " < "];
   string = StringJoin["y", relation, ToString[y[i]]]]]

findrelations[n_] := Module[{},
  areanumber = n;
  onearea = unique2[[areanumber]];
  connections = Partition[First /@ onearea, 2, 1];
  strings = DeleteDuplicates@
    Map[makerelations[#, areanumber] &, connections];
  StringJoin["Area ", ToString[areanumber],
   If[areanumber > 9, ": ", ":  "],
   StringRiffle[strings, " &&\n         "]]]

Show[Plot[Evaluate@Array[y, 4], {x, -1, 1.5},
  PlotLegends -> "Expressions", Axes -> None],
 Graphics[Array[Text[#, centroid[#]] &, len2]]]

Column@Array[findrelations, len2]

Area 1:  y >= 28 - 30 x &&
         y < 3 + 5 x
Area 2:  y >= 2 + 60 x &&
         y >= 28 - 30 x &&
         y < 7 + 90 x
Area 3:  y < 28 - 30 x &&
         y < 7 + 90 x &&
         y < 2 + 60 x &&
         y < 3 + 5 x
Area 4:  y >= 3 + 5 x &&
         y >= 28 - 30 x &&
         y < 2 + 60 x
Area 5:  y < 3 + 5 x &&
         y >= 2 + 60 x &&
         y < 7 + 90 x
Area 6:  y < 28 - 30 x &&
         y >= 2 + 60 x &&
         y >= 3 + 5 x &&
         y < 7 + 90 x
Area 7:  y < 28 - 30 x &&
         y >= 3 + 5 x &&
         y < 2 + 60 x
Area 8:  y < 28 - 30 x &&
         y >= 7 + 90 x &&
         y >= 3 + 5 x
Area 9:  y < 2 + 60 x &&
         y >= 7 + 90 x
Area 10: y >= 28 - 30 x &&
         y >= 7 + 90 x
Area 11: y < 3 + 5 x &&
         y >= 7 + 90 x &&
         y >= 2 + 60 x
qnakjoqk

qnakjoqk3#

完整的Matlab解决方案。120行7021个分区,在0.4秒内完成原始计算(绘图+2秒)。
最后一次。

clear
tic
%rng(15) % fix rng for debug
NLines=9;
DRAW_FINAL_POLY= true==true ; %false true
DRAW_FINALTEXT= true==true ; %false true
DRAW_LINES= false==true; %false true
DRAW_DEBUG= false==true; %false true

% part a - generate lines
NORM_EPS=1e-10;
Lines=10*(rand(NLines,4)-0.5); %(x1,y1,x2,y2)
Color=rand(NLines,3);
Lines(1,:)=[-10,0,+10,0];% x axis if we want to add asix as lines
Lines(2,:)=[0,-10,0,+10];% y axis
Color(1,:)=[1,0,0];% x axis red
Color(2,:)=[0,1,0];% y axis green
Color(3,:)=[0,0,1];% third blue

AllPairs=sortrows(combnk(1:NLines,2));
NPairs=size(AllPairs,1);
[px,py,isok,d] = LineIntersection(Lines(AllPairs(:,1),:) , Lines(AllPairs(:,2),:));
% draw lines and intersections
figure(7); cla; ax=gca;axis(ax,'auto');
if DRAW_LINES
    for iline=1:size(Lines,1)
        line(ax,Lines(iline,[1 3]),Lines(iline,[2 4]),'LineWidth',2','Color',Color(iline,:)); %draw partial line defined by points x1y1x2y2
    end
    for i=1:NPairs %extraploate line to intersection point
        iline1=AllPairs(i,1);
        iline2=AllPairs(i,2);
        line(ax,[Lines(iline1,[1]),px(i)],[Lines(iline1,[2]),py(i)],'LineWidth',2','Color',Color(iline1,:));
        line(ax,[Lines(iline2,[1]),px(i)],[Lines(iline2,[2]),py(i)],'LineWidth',2','Color',Color(iline2,:));
    end
    line(ax,px,py,'LineStyle','none','Marker','*','MarkerSize',8,'Color','g');%draw intersection points
    for i=1:NPairs
        text(px(i)+0.2,py(i)+0.2,sprintf('%d',i),'FontSize',14,'Color','c')
    end
end

% part b - find regions
% 1 for each line sort all its intersections
SortIntr=cell(1,NLines); %(if no parallel lines than NintrsctnsPerline = NLines-1)
IpairiIdx1=cell(1,NLines);
IpairiIdx2=cell(1,NLines);
IpairiIdxI=cell(1,NLines);
for iline=1:NLines
    [idx1,idx2]=find(AllPairs==iline);
    intr=[px(idx1),py(idx1)];
    [~, Isortedintr]=sortrows(intr.*[1,1]); %sort by increaing x y. (prepare clockwise travers)
    SortIntr{iline}=intr(Isortedintr,:);
    IpairiIdx1{iline}=idx1(Isortedintr);
    IpairiIdx2{iline}=3-idx2(Isortedintr); %keep indexes of second line
    IpairiIdxI{iline}=Isortedintr;
end
% 2 traverse from every point and find next closest intersctionn
% we go clockwise x+x-y
PointsInPartition={};
SolIndexInPartition={};
LineIndexInPartition={};
Line4PointsInPartition={};
VisitedSequence=false(NPairs,NPairs); %skip same sequence
count_added=0; count_skipped=0;

for iline=1:NLines
    for ipoint_idx=1:length(SortIntr{iline})-1 %cant start from last point in line
        ipoint=SortIntr{iline}(ipoint_idx,:);
        if DRAW_DEBUG
            delete(findall(ax,'Tag','tmppoint'));
            line(ax,ipoint(1),ipoint(2),'LineStyle','none','Marker','O','MarkerSize',12,'Color','r','tag','tmppoint');%draw intersection points
        end
        current_line=iline;
        isol_idx=IpairiIdx1{current_line}(ipoint_idx);
        current_p_idx=ipoint_idx;
        current_l_next_p_idx=current_p_idx+1;
        next_line=AllPairs(IpairiIdx1{current_line}(current_l_next_p_idx), IpairiIdx2{current_line}(current_l_next_p_idx));
        % next_point_idx = find(IpairiIdx1{current_line}(next_point_idx)==IpairiIdx1{next_line});
        sol_idx_list=[isol_idx]; 
%       if ismember(isol_idx,[7,8,12]),keyboard;end
        point_list=[ipoint];
        line_list=[iline];
        while next_line~=iline
            if DRAW_DEBUG
                delete(findall(ax,'Tag','tmpline'));
                line(ax,Lines(current_line,[1 3]),Lines(current_line,[2 4]),'LineWidth',4','Color',[ 0,0,0 ],'Tag','tmpline');
                line(ax,Lines(next_line,[1 3]),Lines(next_line,[2 4]),'LineWidth',4','Color',[ 1,1,1 ],'Tag','tmpline');
            end
            current_sol_idx=IpairiIdx1{current_line}(current_l_next_p_idx);
            current_p_idx = find(IpairiIdx1{next_line}==current_sol_idx);
            current_line=next_line;
            current_point=SortIntr{current_line}(current_p_idx,:);
            current_nrm=norm(current_point-ipoint);
            current_o=atan2d(-current_point(2)+ipoint(2),current_point(1)-ipoint(1));

            sol_idx_list(end+1)=current_sol_idx; %#ok<SAGROW>
            point_list(end+1,:)=current_point; %#ok<SAGROW>
            line_list(end+1)=current_line;
            if DRAW_DEBUG,line(ax,current_point(1),current_point(2),'LineStyle','none','Marker','O','MarkerSize',12,'Color','m','tag','tmppoint');end %draw intersection points
            %select between two options. next clockwise point is the one with higher angle
            if current_p_idx+1<=length(SortIntr{current_line}) && current_p_idx-1>0
                next_point_1=SortIntr{current_line}(current_p_idx+1,:);
                next_point_2=SortIntr{current_line}(current_p_idx-1,:);
                if norm(next_point_1-ipoint)<NORM_EPS
                    current_l_next_p_idx=current_p_idx+1;
                elseif norm(next_point_2-ipoint)<NORM_EPS
                    current_l_next_p_idx=current_p_idx-1;
                else
                    o1=atan2d(-next_point_1(2)+ipoint(2),next_point_1(1)-ipoint(1));
                    o2=atan2d(-next_point_2(2)+ipoint(2),next_point_2(1)-ipoint(1));
                    if o1>o2
                        current_l_next_p_idx=current_p_idx+1;
                    else
                        current_l_next_p_idx=current_p_idx-1;
                    end
                end
            elseif current_p_idx-1>0
                current_l_next_p_idx=current_p_idx-1;
            else
                current_l_next_p_idx=current_p_idx+1;
            end

            next_p=SortIntr{current_line}(current_l_next_p_idx,:);
            next_o=atan2d(-next_p(2)+ipoint(2),next_p(1)-ipoint(1));
            next_nrm=norm(next_p-ipoint);
            if DRAW_DEBUG,disp([current_nrm,next_nrm,current_o,next_o]);end
            next_line=AllPairs(IpairiIdx1{current_line}(current_l_next_p_idx), IpairiIdx2{current_line}(current_l_next_p_idx));
            next_sol_idx=IpairiIdx1{current_line}(current_l_next_p_idx);
            if VisitedSequence(current_sol_idx,next_sol_idx)
                next_line=-2;
                if DRAW_DEBUG,disp('seq visited');end
                break;
            end
            if next_nrm>NORM_EPS && next_o<current_o || next_o-current_o>180
                next_line=-2;
                if DRAW_DEBUG,disp('next_o<current_o');end
                break; %non clockwise
            end
            assert(next_nrm<NORM_EPS && next_line==iline || next_nrm>=NORM_EPS && next_line~=iline);
        end

        if next_line==iline
            sol_idx_list(end+1)=next_sol_idx; %#ok<SAGROW>
            point_list(end+1,:)=next_p; %#ok<SAGROW>
            PointsInPartition{end+1}=point_list; %#ok<SAGROW>
            SolIndexInPartition{end+1}=sol_idx_list; %#ok<SAGROW>
            %Line4PointsInPartition{end+1}=(Lines(AllPairs(sol_idx_list,1),:), [1; 0])+kron(Lines(AllPairs(sol_idx_list,2),:), [0; 1]);%#ok<SAGROW>
            Line4PointsInPartition{end+1}=Lines(line_list,:);%#ok<SAGROW>

            for i=1:length(sol_idx_list)-1
                VisitedSequence(sol_idx_list(i),sol_idx_list(i+1))=true;
            end
            count_added=count_added+1;
        else
            count_skipped=count_skipped+1;
        end
        if DRAW_DEBUG, disp([next_line==iline, count_added,count_skipped]);end
    end
end

% draw all segments
if DRAW_DEBUG %clear debug
    delete(findall(ax,'Tag','tmppoint'));
    delete(findall(ax,'Tag','tmpline'));
end
NPartition=length(PointsInPartition);
s=sprintf('Lines=%d, Segments=%d, RunTime=%1.2fsec',NLines,NPartition,toc);
title(ax,s);
fprintf([s,newline]);
if DRAW_FINAL_POLY
    hold(ax,'on');
    for i=1:NPartition
        plist=PointsInPartition{i};
        patch(plist(:,1),plist(:,2),i,'FaceAlpha',.3)
        [cx,cy]=ploygon_centroid(plist(:,1),plist(:,2));
        if(DRAW_FINALTEXT),text(cx,cy,sprintf('%d',i),'FontSize',12,'Color','k');end
    end
end

function [px,py,isok,d] = LineIntersection(Line1, Line2)
    %Line1=[x1,y1,x2,y2] Line2=[x3,y3,x4,y4]
    x1=Line1(:,1); y1=Line1(:,2); x2=Line1(:,3); y2=Line1(:,4);
    x3=Line2(:,1); y3=Line2(:,2); x4=Line2(:,3); y4=Line2(:,4);
    d=(x1-x2).*(y3-y4)-(y1-y2).*(x3-x4);%determinant
    px0=(x1.*y2-y1.*x2).*(x3-x4)-(x1-x2).*(x3.*y4-y3.*x4);
    py0=(x1.*y2-y1.*x2).*(y3-y4)-(y1-y2).*(x3.*y4-y3.*x4);
    isok=abs(d)>1e-6;
    px=px0./d;
    py=py0./d;
end

function [xc,yc] = ploygon_centroid(x,y)
    xs = circshift(x,-1);
    ys = circshift(y,-1);
    area = 0.5*sum (x.*ys-xs.*y);
    xc = sum((x.*ys-xs.*y).*(x+xs))/(6*area);
    yc = sum((x.*ys-xs.*y).*(y+ys))/(6*area);
end

第一次

bejyjqdl

bejyjqdl4#

EDIT New solution, using itertools. Old solution below. Using the python package itertools , I believe this solution is faster than the original. However, it is still extremely slow and is not feasible beyond about 20-ish lines. For 120 lines, none of these will terminate.

import itertools

strings = ["<=", ">="]

fxs = ["5x + 3", "-60x + 7", ...]

parts = []

if __name__ == "__main__":
    print(len(fxs))
    out = itertools.product(strings, repeat=len(fxs))
    for i in out:
        curpart = ""
        for j in range(len(i)):
            print(curpart)
            if j != len(i):
                curpart = curpart + "y " + i[j] + fxs[j] + " && "
            else:
                curpart = curpart + "y " + i[j] + fxs[j]
        parts.append(curpart)
    print(parts)

OLD SOLUTION

Here's what ended up working for me, in python:

fxs = [list of functions as strings, ex: "5x+3","6x+4"... etc.]

parts = []

if __name__ == "__main__":
    for f in fxs:
        fhalf = "y >= " + f
        shalf = "y <= " + f
        if len(parts) == 0:
            parts.append(fhalf)
            parts.append(shalf)
        else:
            parts1 = [s + " && " + fhalf for s in parts]
            parts2 = [s + " && " + shalf for s in parts]
            parts = parts1 + parts2
    print(parts)

For the above example, this outputs:

['y >= 5x+3 && y >= 90x+7 && y >= -30x+28 && y >= 60x+2', 'y <= 5x+3 && y >= 90x+7 && y >= -30x+28 && y >= 60x+2', 'y >= 5x+3 && y <= 90x+7 && y >= -30x+28 && y >= 60x+2', 'y <= 5x+3 && y <= 90x+7 && y >= -30x+28 && y >= 60x+2', 'y >= 5x+3 && y >= 90x+7 && y <= -30x+28 && y >= 60x+2', 'y <= 5x+3 && y >= 90x+7 && y <= -30x+28 && y >= 60x+2', 'y >= 5x+3 && y <= 90x+7 && y <= -30x+28 && y >= 60x+2', 'y <= 5x+3 && y <= 90x+7 && y <= -30x+28 && y >= 60x+2', 'y >= 5x+3 && y >= 90x+7 && y >= -30x+28 && y <= 60x+2', 'y <= 5x+3 && y >= 90x+7 && y >= -30x+28 && y <= 60x+2', 'y >= 5x+3 && y <= 90x+7 && y >= -30x+28 && y <= 60x+2', 'y <= 5x+3 && y <= 90x+7 && y >= -30x+28 && y <= 60x+2', 'y >= 5x+3 && y >= 90x+7 && y <= -30x+28 && y <= 60x+2', 'y <= 5x+3 && y >= 90x+7 && y <= -30x+28 && y <= 60x+2', 'y >= 5x+3 && y <= 90x+7 && y <= -30x+28 && y <= 60x+2', 'y <= 5x+3 && y <= 90x+7 && y <= -30x+28 && y <= 60x+2']

The code is extremely simple and the output is exactly as desired.
The idea is basically that each line divides the plane into two halves, and so the loop (which iterates over the set of lines) intersects each partition already found with each half and adds the two new partitions to the set, while removing the original, in-intersected partition. The result doesn't always give the simplest possible condition (some of the conditions may be redundant), but they do give a complete description of each partition.
While this solution works, with 120 lines, it is pretty slow. I'd be interested in seeing if there are more efficient ways to accomplish this, using this method or otherwise.

mxg2im7a

mxg2im7a5#

到目前为止所建议的方法的另一种方法是使用三角形。
从一组三角形开始(可能包含无穷远的点),它们的并集是一个覆盖整个平面的多边形。例如,您可以开始使用作为算法输入的前两条直线分割平面所得到的四个三角形,或者使用由原点O、OX轴的正边和240度并且也从O开始的两条半线。
现在,对于每一条作为输入的线,您都要查找它所穿过的三角形,对于每一个三角形,您都要在与线的交点处将其截断,并将其替换为三个三角形(或者在退化的情况下是两个)覆盖相同区域的新三角形。然后你必须找到那些旧的相交三角形是其一部分的多边形,并且用两个多边形来替换其中的每一个,这两个多边形是根据它们所在的线的边来分割它所包含的三角形而得到的。
使用三角形的好处是它使计算更容易,你可以把它放到一些索引数据结构中,比如k-d树,这样可以更快地计算哪些多边形被某条线切割,你可以用一个复杂度为O(NlogN)的算法来结束,其中N是形成平面划分的多边形的数量。

相关问题