Use segment tree and scan line algorithm to calculate polygon ROI when self-intersection occurs

[IOSurfer]

I implement my own version of it to calculate area,mean and standard deviation of the ROI.It performs quite well in real-time.

export class Aabb2D {
  xMin: number;
  yMin: number;
  xMax: number;
  yMax: number;

  constructor(x0: number, y0: number, x1: number, y1: number) {
    this.xMin = Math.min(x0, x1);
    this.xMax = Math.max(x0, x1);
    this.yMin = Math.min(y0, y1);
    this.yMax = Math.max(y0, y1);
  }

  Union(bound: Aabb2D): Aabb2D {
    return new Aabb2D(
      Math.min(this.xMin, bound.xMin),
      Math.min(this.yMin, bound.yMin),
      Math.max(this.xMax, bound.xMax),
      Math.max(this.yMax, bound.yMax)
    );
  }
}

import { Aabb2D } from "./Aabb2D";

export class Segment {
  x0: number;
  y0: number;
  x1: number;
  y1: number;
  xMin: number;
  xMax: number;

  constructor(x0: number, y0: number, x1: number, y1: number) {
    this.x0 = x0;
    this.y0 = y0;
    this.x1 = x1;
    this.y1 = y1;
    this.xMin = Math.min(x0, x1);
    this.xMax = Math.max(x0, x1);
  }

  GetBound(): Aabb2D {
    return new Aabb2D(this.x0, this.y0, this.x1, this.y1);
  }
}

export class SegmentTreeNode {
  leftChildIndex: number;
  rightChildIndex: number;
  start: number;
  end: number;
  isLeafNode: boolean;
  segmentIndices: number[];

  constructor(start: number = 0, end: number = 0) {
    this.leftChildIndex = -1;
    this.rightChildIndex = -1;
    this.start = start;
    this.end = end;
    this.isLeafNode = false;
    this.segmentIndices = [];
  }
}


import { Segment } from './Segment';
import { SegmentTreeNode } from './SegmentTreeNode';
import { Aabb2D } from './Aabb2D';

export class SegmentTree {
  segmentList_: Segment[];
  nodeList_: SegmentTreeNode[];
  start_: number;
  end_: number;
  length_: number;
  rootIndex_: number;

  constructor(segmentList_: Segment[]) {
    this.segmentList_ = segmentList_;
    this.nodeList_ = [];
    this.start_ = 0;
    this.end_ = 0;
    this.length_ = 0;
    this.rootIndex_ = -1;
  }

  buildTree() {
    const vertices: number[] = [];

    this.segmentList_.forEach(({ x0 }) => vertices.push(x0));

    vertices.sort((a, b) => a - b);

    const mergeThreshold = 5.0;

    const newVertices: number[] = [vertices[0]];
    for (let i = 1; i < vertices.length - 1; ++i) {
      if (
        Math.abs(vertices[i] - newVertices[newVertices.length - 1]) >
        mergeThreshold
      ) {
        newVertices.push(vertices[i]);
      }
    }
    newVertices.push(vertices[vertices.length - 1]);

    this.start_ = newVertices[0];
    this.end_ = newVertices[newVertices.length - 1];
    this.length_ = this.end_ - this.start_;

    let nodeCnt = 0;
    const leafNum = newVertices.length - 1;

    for (let i = 1; i < newVertices.length; i++) {
      const node: SegmentTreeNode = {
        leftChildIndex: -1,
        rightChildIndex: -1,
        start: newVertices[i - 1],
        end: newVertices[i],
        isLeafNode: true,
        segmentIndices: [],
      };
      this.nodeList_.push(node);
      nodeCnt++;
    }

    // Insert segments non-recursively
    for (
      let curSegIndex = 0;
      curSegIndex < this.segmentList_.length;
      curSegIndex++
    ) {
      let curNodeIndex = 0;
      while (
        curNodeIndex < leafNum &&
        this.nodeList_[curNodeIndex].end <= this.segmentList_[curSegIndex].xMin
      ) {
        curNodeIndex++;
      }
      while (
        curNodeIndex < leafNum &&
        this.nodeList_[curNodeIndex].start <=
          this.segmentList_[curSegIndex].xMax
      ) {
        this.nodeList_[curNodeIndex].segmentIndices.push(curSegIndex);
        curNodeIndex++;
      }
    }

    this.rootIndex_ = this.mergeSegment(nodeCnt, 0, leafNum - 1);
  }

  mergeSegment(nodeCnt: number, startIndex: number, endIndex: number): number {
    const stack: {
      startIndex: number;
      endIndex: number;
      returning: boolean;
    }[] = [];
    const results: number[] = [];

    stack.push({ startIndex, endIndex, returning: false });

    while (stack.length > 0) {
      const cur = stack.pop()!;
      if (cur.startIndex >= cur.endIndex) {
        results.push(cur.endIndex);
        continue;
      }
      if (!cur.returning) {
        stack.push({
          startIndex: cur.startIndex,
          endIndex: cur.endIndex,
          returning: true,
        });
        const midIndex = Math.floor((cur.startIndex + cur.endIndex) / 2);
        stack.push({
          startIndex: cur.startIndex,
          endIndex: midIndex,
          returning: false,
        });
        stack.push({
          startIndex: midIndex + 1,
          endIndex: cur.endIndex,
          returning: false,
        });
      } else {
        const leftIndex = results.pop()!;
        const rightIndex = results.pop()!;
        const node: SegmentTreeNode = {
          start: this.nodeList_[leftIndex].start,
          end: this.nodeList_[rightIndex].end,
          isLeafNode: false,
          leftChildIndex: leftIndex,
          rightChildIndex: rightIndex,
          segmentIndices: [],
        };
        this.nodeList_.push(node);
        results.push(nodeCnt);
        nodeCnt++;
      }
    }

    return nodeCnt - 1;
  }

  query(x: number, results: number[]): number {
    results.length = 0;
    let curNodeIndex = this.rootIndex_;

    while (!this.nodeList_[curNodeIndex].isLeafNode) {
      if (
        this.nodeList_[curNodeIndex].start > x ||
        this.nodeList_[curNodeIndex].end <= x
      ) {
        return 0;
      }

      const { leftChildIndex } = this.nodeList_[curNodeIndex];
      const { rightChildIndex } = this.nodeList_[curNodeIndex];

      if (
        this.nodeList_[leftChildIndex].start > x ||
        this.nodeList_[leftChildIndex].end <= x
      ) {
        curNodeIndex = rightChildIndex;
        continue;
      } else {
        curNodeIndex = leftChildIndex;
      }
    }

    if (
      this.nodeList_[curNodeIndex].start <= x &&
      this.nodeList_[curNodeIndex].end > x
    ) {
      this.nodeList_[curNodeIndex].segmentIndices.forEach((itr) => {
        const seg = this.segmentList_[itr];
        if (x >= seg.xMin && x < seg.xMax) {
          results.push(
            seg.y0 + ((seg.y1 - seg.y0) / (seg.x1 - seg.x0)) * (x - seg.x0)
          );
        }
      });

      results.sort((a, b) => a - b);
      if (results.length % 2 !== 0) {
        results.pop();
      }

      let deltaY = 0;
      for (let i = 1; i < results.length; i += 2) {
        deltaY += results[i] - results[i - 1];
      }

      return deltaY;
    }
    return 0;
  }

  queryRange(start: number, end: number, results: number[][]): number {
    let deltaY = 0;
    for (let x = start; x <= end; x++) {
      const result: number[] = [];
      results.push(result);
      deltaY += this.query(x, results[x - start]);
    }
    return deltaY;
  }

  querySamplesRange(
    sampleNum: number,
    start: number,
    end: number,
    results: number[][]
  ): number {
    let deltaY = 0;
    const stride = (end - start) / (sampleNum - 1.0);
    for (let i = 0; i < sampleNum; i++) {
      const result: number[] = [];
      results.push(result);
      deltaY += stride * this.query(start + i * stride, results[i]);
    }
    return deltaY;
  }

  print(): void {
    const q: number[] = [this.rootIndex_];

    while (q.length) {
      const index = q.shift()!;
      const node = this.nodeList_[index];

      console.log(`Index: ${index}`);
      console.log(`Start: ${node.start}`);
      console.log(`End: ${node.end}`);
      console.log(`Left: ${node.leftChildIndex}`);
      console.log(`Right: ${node.rightChildIndex}`);
      console.log(`Leaf: ${node.isLeafNode}`);

      console.log('Segment:', node.segmentIndices.join(';'));

      if (node.leftChildIndex !== -1) {
        q.push(node.leftChildIndex);
      }
      if (node.rightChildIndex !== -1) {
        q.push(node.rightChildIndex);
      }

      console.log('\n');
    }
  }
}

export function printSegments(segments: Segment[]) {
  segments.forEach((seg, i) => {
    console.log(`Segment: ${i} (${seg.x0}, ${seg.y0}) (${seg.x1}, ${seg.y1})`);
  });
}

export function printResults(
  area: number,
  results: number[][],
  verbose = false
) {
  console.log(`Area: ${area}`);
  if (verbose) {
    results.forEach((result, i) => {
      console.log(`Scan ${i}: ${result.join('; ')}`);
    });
  }
}

export function getSegmentsBound(segments: Segment[]): Aabb2D {
  if (segments.length === 0) {
    throw new Error('Segment list cannot be empty.');
  }

  let bound = new Aabb2D(
    segments[0].x0,
    segments[0].y0,
    segments[0].x1,
    segments[0].y1
  );

  segments.forEach((seg) => {
    const segBound = new Aabb2D(seg.x0, seg.y0, seg.x1, seg.y1);
    bound = bound.Union(segBound);
  });

  return bound;
}

function calculate(id: PathID) {
    if (!(id in pathByID.value)) {
      console.warn(`Path with ID ${id} does not exist.`);
      return;
    }

    const path = pathByID.value[id];

    if (!path.pointsInIndex) {
      console.warn(`Path with ID ${id} does not have a pointsInIndex array.`);
      return;
    }

    const points = path.pointsInIndex;

    if (path.placed && points && points.length >= 3) {
      const segments: Segment[] = [];

      for (let x = 1; x < points.length; x++) {
        segments.push(
          new Segment(
            points[x - 1][0],
            points[x - 1][1],
            points[x][0],
            points[x][1]
          )
        );
      }
      segments.push(
        new Segment(
          points[points.length - 1][0],
          points[points.length - 1][1],
          points[0][0],
          points[0][1]
        )
      );

      const segmentTree = new SegmentTree(segments);
      segmentTree.buildTree();

      const bound: Aabb2D = getSegmentsBound(segments);
      const queryResults: number[][] = [];

      const minSampleNum = 50;
      const sampleNum = Math.max(
        Math.ceil(bound.xMax) - Math.floor(bound.xMin) + 1,
        minSampleNum
      );
      const stride =
        (Math.ceil(bound.xMax) - Math.floor(bound.xMin)) / (sampleNum - 1.0);

      const area = segmentTree.querySamplesRange(
        sampleNum,
        Math.floor(bound.xMin),
        Math.ceil(bound.xMax),
        queryResults
      );

      const imageData = imageSotre.dataIndex[path.imageID];
      const extent = imageData.getExtent();
      // Welford’s Algorithm
      let sampleCnt = 0;
      let mean = 0;
      let m2 = 0;
      for (let i = 0; i < queryResults.length; i++) {
        const curX = Math.floor(bound.xMin) + i * stride;
        for (let j = 1; j < queryResults[i].length; j += 2) {
          let curY = queryResults[i][j - 1];
          while (curY <= queryResults[i][j]) {
            curY += stride;
            const index = [];
            for (let component = 0, flag = 0; component < 3; component++) {
              if (component === path.lookAxis) {
                index.push(path.slice);
              } else if (flag === 0) {
                index.push(Math.round(curX));
                flag = 1;
              } else {
                index.push(Math.round(curY));
              }
            }
            let offsetIndex = imageData.computeOffsetIndex([
              index[0],
              index[1],
              index[2],
            ]);
            let result = Number.NaN;
            for (let component = 0; component < 3; component++) {
              if (
                index[component] < extent[2 * component] ||
                index[component] > extent[2 * component + 1]
              ) {
                offsetIndex = Number.NaN;
              }
            }
            if (!Number.isNaN(offsetIndex)) {
              result = imageData
                .getPointData()
                .getScalars()
                .getComponent(offsetIndex, 0);
              sampleCnt++;
              const delta = result - mean;
              mean += delta / sampleCnt;
              m2 += delta * (result - mean);
            }
          }
        }
      }

      mean = sampleCnt > 1 ? mean : 0;
      // let $E(\hat{\sigma}) = \sigma$
      const stdDev = sampleCnt > 1 ? Math.sqrt(m2 / (sampleCnt - 1)) : 0;
      console.log('area:', area, 'mean', mean, 'stdDev', stdDev);
      updatePath(id, { areaInPixel: area, mean, stdDev });
    }
  }

[IOSurfer]

I initially implemented a C++ version of it, but I found it works quite fast in typescript, so I did not compile it into a WASM module.

[PaulHax]

This is very impressive. Thank you for sharing.