ถอดแบบจาก PDF ที่วาดไม่ตรงสเกล แล้ววาดกลับให้ตรงระยะ

พอดีผมได้แบบผังโรงงานมาเป็น PDF ไฟล์เดียว โจทย์คือขึ้นโมเดลสามมิติให้ตรงระยะจริง ทั้งตัวห้องและของที่วางอยู่ข้างใน ตอนแรกคิดว่าง่าย แค่วัด px แล้วคูณสเกล จบ ปรากฏว่าพอทำจริงแล้วของทุกชิ้นเพี้ยน แท่งที่ควรจะเป็นก้อนหนากลายเป็นแผ่นบาง วงกลมกลายเป็นวงรี

สาเหตุไม่ได้อยู่ที่โค้ด แต่อยู่ที่ตัวแบบเอง ผมเลยเขียนบันทึกวิธีทำทั้งชุดเก็บไว้ ตั้งแต่การหาสเกลจริงของแต่ละห้อง ไปจนถึงการนับจำนวนชิ้นแบบที่เถียงไม่ได้ และการเอาตัวเลขวาดกลับลงโมเดล

เรื่องดึงขอบเขตแปลงที่ดินจากระบบแผนที่ราชการ ผมเขียนแยกไว้แล้วที่ นำเข้าขอบเขตแปลงที่ดินเข้า SketchUp อันนั้นคือข้อมูลที่เป็นตัวเลขอยู่แล้ว ส่วนอันนี้คืองานอีกแบบหนึ่ง คือมีแต่ภาพกระดาษ ต้องคั้นตัวเลขออกมาเอง

แบบที่วาดสวย ไม่ได้แปลว่าวัดได้

แบบแนวนำเสนอ (presentation) กับแบบก่อสร้างจริงคนละเรื่องกัน แบบนำเสนอเขาวาดให้ "เล่าเรื่องรู้เรื่อง" ห้องที่มีของเยอะเขาก็ยืดให้ใหญ่ขึ้นเพื่อให้ใส่สัญลักษณ์ได้ครบ ห้องที่ว่าง ๆ ก็บีบให้แคบลงเพื่อประหยัดกระดาษ ระยะที่เขียนกำกับไว้เป็นระยะจริง แต่รูปที่วาดไม่ได้ยึดตามระยะนั้น

ผมวัดสเกลจริงของแต่ละห้องออกมาจากไฟล์เดียวกัน ที่ 300 dpi หน่วยเป็น pixel ต่อเมตร ตัวเลขที่ได้เป็นแบบนี้

ดูให้ดีสองชั้นครับ ชั้นแรกคือในห้องเดียวกัน แกน x กับ y ยังคนละสเกล ห้อง A ต่างกันเกือบเท่าตัว (49.5 กับ 29.8) แปลว่าอะไรที่วาดเป็นวงกลมในห้องนั้น พอวัด px จะได้เป็นวงรีเสมอ ชั้นที่สองคือข้ามห้องแล้วยิ่งไม่เท่ากันไปอีก แกน y ของห้อง B อยู่ที่ 22.0 แต่ห้อง C อยู่ที่ 41.0 ต่างกันเกือบสองเท่า

รูปทรงบนกระดาษหลอกตา ตัวอย่างที่ชัดที่สุด อยู่ในห้อง B ครับ

ในนั้นมีตะแกรงร่อนอยู่ตัวหนึ่ง บนกระดาษมันเป็นแท่งผอมยาว วัด px ได้สัดส่วนประมาณ 4.9 : 1 ถ้าวาดตามที่เห็นก็จะได้แผ่นบาง ๆ ยาว ๆ วางอยู่กลางห้อง

พอเอาสเกลจริงของห้องนั้นมาหาร (แกน x 49.4 แกน y 22.0) ผลออกมาคือ 3.43 x 2.00 เมตร คือสัดส่วน 1.7 : 1 เป็นก้อนหนา ๆ เกือบจัตุรัส คนละอย่างกับที่ตาเห็นเลย

ส่วนตัวผมถือว่านี่คือเหตุผลที่ต้องทำทั้งกระบวนการนี้ ไม่ใช่แค่ "วาดตามรูป" เพราะถ้าวาดตามรูป ของชิ้นนี้ชิ้นเดียวก็ผิดไปแล้วสองเท่ากว่า

ขั้นแรก ซูมดูก่อนจะเชื่ออะไร

ผมใช้ PyMuPDF (fitz) ครอปเฉพาะกรอบห้องที่สนใจแล้วเรนเดอร์ใหม่ที่ความละเอียดสูง ตรงนี้สำคัญมาก เพราะเรนเดอร์ทั้งหน้าที่ 300 dpi แล้วซูมดูในโปรแกรมดูรูป ภาพจะแตกจนนับของไม่ได้ ต้องให้ fitz เรนเดอร์เฉพาะกรอบนั้นที่ dpi สูงตั้งแต่แรก

import fitz

PDF = 'plan/factory_plan.pdf'
DPI = 300                                   # ความละเอียดที่ใช้อ่านพิกัดทั้งไฟล์

doc = fitz.open(PDF)
page = doc[0]

# กรอบที่วัดมาเป็น px ที่ 300 dpi -> ต้องแปลงกลับเป็น point (72 dpi) ให้ fitz
x0, y0, x1, y1 = 1993, 1698, 3072, 2122
r = fitz.Rect(x0 / DPI * 72, y0 / DPI * 72, x1 / DPI * 72, y1 / DPI * 72)

page.get_pixmap(dpi=800, clip=r).save('zoom_room_a.png')

dpi=800 ไม่ได้ใส่มาเล่น ๆ ตอนที่ต้องนับจำนวนชิ้นที่วางชิดกัน ผมต้องขึ้นไปถึงช่วง 600-900 dpi ถึงจะแยกออกว่าอันไหนเป็นชิ้นเดียวกัน อันไหนเป็นสองชิ้นวางติดกัน ต่ำกว่านั้นเส้นขอบจะเชื่อมกันจนนับผิด

หาผนังห้อง เพื่อคาลิเบรตห้องนั้นโดยเฉพาะ

ผนังในแบบมีสองแบบ ผนังทึบวาดเป็นเส้นเต็ม ส่วนขอบเขตพื้นที่บางที่วาดเป็นเส้นประ สองแบบนี้จับด้วยเกณฑ์เดียวกันไม่ได้

ฟังก์ชันที่ผมใช้หาเส้นเป็นแบบนี้ คืนเลขแถวกับเลขคอลัมน์ที่เข้าเกณฑ์ออกมาให้ดูด้วยตาก่อน แล้วค่อยเลือกเอง ไม่ปล่อยให้มันเดาให้

import numpy as np

pix = doc[0].get_pixmap(dpi=DPI)
img = (np.frombuffer(pix.samples, dtype=np.uint8)
         .reshape(pix.height, pix.width, pix.n)[:, :, :3].astype(int))
black = (img.sum(axis=2) < 200)


def walls(y0, y1, x0, x1, frac=0.25):
    """frac ต่ำ = จับเส้นประได้  ·  frac สูง = เอาเฉพาะผนังทึบ"""
    s = black[y0:y1, x0:x1]
    rows = [i + y0 for i, v in enumerate(s.sum(axis=1)) if v > (x1 - x0) * frac]
    cols = [i + x0 for i, v in enumerate(s.sum(axis=0)) if v > (y1 - y0) * frac]
    return rows, cols

พอได้ขอบซ้าย-ขวา-บน-ล่างของห้อง (L, R, T, B เป็น px) แล้ว ผมผูกมันกับระยะจริงที่อ่านจากโมเดล ไม่ใช่ตัวเลขที่เขียนกำกับบนแบบ เพราะแบบเขียนระยะรวมผิดได้บ่อยกว่าที่คิด ส่วนโมเดลเป็นของที่เราตรวจแล้วว่าถูก

ROOM  = dict(L=1993.5, R=3072.0, T=1698.0, B=2122.0)   # ขอบผนังห้อง (px @300dpi)
MODEL = dict(x0=75.0, x1=87.0, y0=100.0, y1=107.0)     # ระยะจริงจากโมเดล (เมตร)

L, R, T, B = ROOM['L'], ROOM['R'], ROOM['T'], ROOM['B']
W = MODEL['x1'] - MODEL['x0']
H = MODEL['y1'] - MODEL['y0']

SX, SY = (R - L) / W, (B - T) / H       # px ต่อเมตร คนละค่ากันได้ = แบบบิด

X = lambda px: round(MODEL['x0'] + (px - L) / (R - L) * W, 2)
Y = lambda py: round(MODEL['y1'] - (py - T) / (B - T) * H, 2)

สองบรรทัดล่างคือหัวใจทั้งหมดของงานนี้ครับ X กับ Y แปลง px เป็นเมตรของห้องนั้น ๆ และเปลี่ยนค่าตามห้อง แกน y กลับด้านเพราะ pixel นับจากบนลงล่าง แต่พิกัดโมเดลนับจากล่างขึ้นบน

นับจำนวนชิ้นด้วย colour-blob ไม่ใช่ด้วยตา

แบบชุดนี้ระบายสีของแต่ละชนิดไว้แล้ว ถังสีเขียว รางสีฟ้า สายพานสีส้ม ตู้สีเทา ตรงนี้เป็นของขวัญเลย เพราะแยกชนิดได้จากสีตรง ๆ แต่ห้ามเดาค่าสี ให้ดึง histogram สีจริงในกรอบห้องนั้นมาดูก่อนเสมอ สีที่เห็นบนจอกับค่าที่อยู่ในไฟล์ไม่ตรงกันบ่อยมาก โดยเฉพาะถ้า PDF ผ่านการบีบมา

from collections import Counter, deque


def palette(y0, y1, x0, x1, n=10):
    """ดูสีจริงในกรอบห้องก่อน อย่าเดาค่าสี"""
    sub = img[y0:y1, x0:x1].reshape(-1, 3)
    return Counter(map(tuple, sub)).most_common(n)


def mask(rgb, tol=60):
    return np.abs(img - np.array(rgb)).sum(axis=2) < tol

ได้ mask ของสีแล้วก็ BFS หา connected component ธรรมดา ๆ ตัวนี้แหละที่ทำให้ "จำนวนชิ้น" กลายเป็นตัวเลขที่เถียงไม่ได้ แทนที่จะเป็นการเพ่งจอแล้วนับ

def blobs(m, y0, y1, x0, x1, minpx=100):
    sub = m[y0:y1, x0:x1]
    h, w = sub.shape
    vis = np.zeros((h, w), bool)
    out = []
    for i in range(h):
        for j in range(w):
            if not sub[i, j] or vis[i, j]:
                continue
            q = deque([(i, j)]); vis[i, j] = True; pts = []
            while q:
                a, b = q.popleft(); pts.append((a, b))
                for da, db in ((1, 0), (-1, 0), (0, 1), (0, -1)):
                    na, nb = a + da, b + db
                    if 0 <= na < h and 0 <= nb < w and sub[na, nb] and not vis[na, nb]:
                        vis[na, nb] = True
                        q.append((na, nb))
            if len(pts) >= minpx:
                ar = np.array(pts)
                out.append(dict(n=len(pts),
                                y0=int(ar[:, 0].min()) + y0, y1=int(ar[:, 0].max()) + y0,
                                x0=int(ar[:, 1].min()) + x0, x1=int(ar[:, 1].max()) + x0))
    return out

วงกลมที่ไม่กลม และ blob ที่จริง ๆ เป็นของสองชิ้น

พอได้กรอบของแต่ละ blob มาแล้ว มีสองกับดักที่เจอทุกงาน

กับดักที่หนึ่ง วงกลมที่วัดได้ w ไม่เท่ากับ h อย่าเพิ่งคิดว่าเป็นวงรีจริง ๆ ในแบบ เพราะสเกลแกน x กับ y ไม่เท่ากันอยู่แล้ว วงกลมจริงจะออกมาเป็นวงรีเสมอ วิธีคิดเส้นผ่าศูนย์กลางที่ถูกคือแปลง w และ h เป็นเมตรก่อน แล้วใช้ค่าเฉลี่ยเรขาคณิต sqrt(w * h)

กับดักที่สอง blob ที่สัดส่วนผิดปกติ ถ้าเจอ blob ที่สูงหรือกว้างเป็นสองเท่าของเพื่อนในกลุ่มเดียวกัน แทบทุกครั้งคือของสองชิ้นที่วางชิดกันจน pixel เชื่อมถึงกัน BFS มองเป็นก้อนเดียว ต้องกลับไปซูมดูแล้วแยกเอง

def report(bs, label):
    print('==', label, len(bs))
    for b in bs:
        w = (b['x1'] - b['x0']) / SX          # px -> เมตร ด้วยสเกลแกน x ของห้องนี้
        h = (b['y1'] - b['y0']) / SY          # px -> เมตร ด้วยสเกลแกน y ของห้องนี้
        cx = (X(b['x0']) + X(b['x1'])) / 2.0
        cy = (Y(b['y0']) + Y(b['y1'])) / 2.0
        # วงกลม: สเกลบิดทำให้ w != h เสมอ -> เส้นผ่าศูนย์กลางใช้ค่าเฉลี่ยเรขาคณิต
        print('   c=(%.2f,%.2f)  w=%.2f h=%.2f  dia~%.2f' % (cx, cy, w, h, (w * h) ** 0.5))

เขียนกลับลงโมเดล: สองกับดักของ SketchUp

ได้พิกัดเมตรครบแล้วก็เขียนเป็นสคริปต์ Ruby ยิงเข้าโมเดล ตรงนี้เจอสองเรื่องที่กินเวลาผมไปพอสมควร

เรื่องแรก กลุ่มในโมเดลมี transformation ไม่เป็น identity ห้องที่ลูกค้าเคยขยับมาก่อน จะมี origin ของกลุ่มเลื่อนไปแล้ว ถ้าใส่พิกัดโลกลงไปตรง ๆ ของจะเลื่อนตามไปหมด ต้องแปลงเป็นพิกัด local ของกลุ่มนั้นทุกจุด

# ส่งพิกัดโลก (เมตร) เข้ามาเสมอ แล้วให้ lpt แปลงเป็น local ของกลุ่มปลายทางให้
def lpt(parent, x, y, z = 0)
  Geom::Point3d.new(x.m, y.m, z.m).transform(parent.transformation.inverse)
end

def box(parent, x0, y0, x1, y1, rgb, name, z = 0)
  g = parent.entities.add_group
  g.name = name
  f = g.entities.add_face([lpt(parent, x0, y0, z), lpt(parent, x1, y0, z),
                           lpt(parent, x1, y1, z), lpt(parent, x0, y1, z)])
  f.reverse! if f.normal.z < 0
  if rgb
    f.material = rgb
    f.back_material = rgb
  end
  g
end

เช็คก่อนเริ่มทุกครั้งด้วย g.transformation.identity? ถ้าคืน false แปลว่าต้องใช้ lpt แน่นอน

เรื่องที่สอง ทุกอย่างอยู่ที่ z ศูนย์ระนาบเดียวกัน ผลคือ z-fighting ภาพกะพริบตอนหมุนมุมมอง และมองไม่ออกว่าอะไรทับอะไร ผมแก้ด้วยการยกชิ้นที่ในแบบอยู่ "ข้างบน" ขึ้น 2 เซนติเมตร และกดชิ้นที่ถูกบังลง 2 เซนติเมตร แค่นี้ภาพก็นิ่งและอ่านลำดับชั้นได้ตรงกับแบบ

1. เขียนสคริปต์ให้ idempotent คือลบกลุ่มย่อยที่ชื่อของตัวเองทิ้งก่อน แล้วค่อยสร้างใหม่ รันซ้ำกี่รอบก็ได้ผลเดิม ไม่ทับซ้อน
2. บรรทัดสุดท้ายของสคริปต์คืน status string ที่บอกจำนวนชิ้นแต่ละชนิดออกมา เอาไว้เทียบกับจำนวนที่ BFS นับได้
3. อย่าใส่ Sketchup::Text ลงในโมเดล มันพังตอน zoom_extents ป้ายหมายเลขให้วาดเป็นวงกลมเส้นโปร่งแทน (สร้าง face แล้วลบ face ทิ้ง เหลือแต่เส้น)
4. ยกหรือกด z ของชิ้นที่ทับกันตามลำดับชั้นที่เห็นในแบบ

ตรวจแบบที่เถียงไม่ได้: warp กลับเป็นสเกลจริงแล้ววางทับ

ตรวจด้วยการมองรูปเทียบรูปไม่ได้ผลครับ เพราะแบบมันบิดอยู่แล้ว ภาพจากโมเดลกับภาพจากกระดาษหน้าตาต่างกันเป็นธรรมดา ต่อให้โมเดลถูกก็ยังดูไม่เหมือน

วิธีที่ใช้ได้จริงคือทำกลับด้าน คือ resample แบบ PDF ทีละบล็อกด้วย px ต่อเมตรของบล็อกนั้น ลงบนผืนผ้าใบที่เป็นหน่วยเมตรร่วมกัน เท่ากับ "ดัดแบบให้ตรงสเกล" แล้วค่อยเอา render จาก SketchUp กรอบโลกเดียวกันมาวางทับ ถ้าตรงแปลว่าโมเดลถูก

from PIL import Image

PPM = 60.0                                   # px ต่อเมตรของภาพผลลัพธ์
X0, X1, Y0, Y1 = 81.0, 97.0, 129.0, 138.0    # กรอบโลก (เมตร) ที่จะ warp
# (Ymin, Ymax, pxAt_Xa, Xa, pxPerM_x, pyAt_Yb, Yb, pyPerM_y) -> หนึ่งบล็อก = หนึ่งสเกล
BANDS = [
    (137.0, 145.0, 2328, 81.0, (3070 - 2328) / 15.0, 619, 145.0, (857 - 619) / 8.0),
    (132.0, 137.0, 2428, 83.0, (3070 - 2428) / 13.0, 887, 137.0, (997 - 887) / 5.0),
]

W, H = int((X1 - X0) * PPM), int((Y1 - Y0) * PPM)
out = np.full((H, W, 3), 255, np.uint8)
for j in range(H):
    Ym = Y1 - (j + .5) / PPM
    b = next((b for b in BANDS if b[0] - .001 <= Ym <= b[1] + .001), None)
    if not b:
        continue
    _, _, pxa, Xa, sx, pyb, Yb, sy = b
    py = int(round(pyb + (Yb - Ym) * sy))
    if not (0 <= py < img.shape[0]):
        continue
    for i in range(W):
        px = int(round(pxa + (X0 + (i + .5) / PPM - Xa) * sx))
        if 0 <= px < img.shape[1]:
            out[j, i] = img[py, px]
Image.fromarray(out).save('warp_true.png')

ฝั่ง SketchUp ก็ตั้งกล้อง ortho ให้กรอบตรงกับ X0-X1 Y0-Y1 ข้างบนเป๊ะ ๆ แล้วเขียนภาพออกมา

v = Sketchup.active_model.active_view
ctr    = Geom::Point3d.new(89.0.m, 133.5.m, 0)          # กลางกรอบเดียวกับที่ warp
ctrtop = Geom::Point3d.new(89.0.m, 133.5.m, 500.m)
v.camera = Sketchup::Camera.new(ctrtop, ctr, Geom::Vector3d.new(0, 1, 0))
v.camera.perspective = false
v.camera.height = 9.0.m                                  # = Y1 - Y0
v.write_image('cmp_model.png', 960, 540, true, 0.0)

พอวางสองภาพซ้อนกันแล้วเส้นตรง ก็จบครับ ถ้าไม่ตรง จะเห็นทันทีว่าไม่ตรงตรงไหนและกี่เมตร ไม่ใช่ความรู้สึกว่า "ดูแปลก ๆ"

กฎที่ผมยึดทุกครั้ง

อีกข้อที่ผมเพิ่มมาทีหลังคือรายงานตัวเลขทุกครั้ง บอกให้ครบว่าอ่านค่ามาจาก px ไหน สเกลเท่าไร ได้พิกัดอะไร และอันไหนที่ "อนุมาน" ไม่ใช่ "วัดมา" เพราะเวลาผลไม่ตรงกับที่ลูกค้าคาด จะได้ไล่กลับได้ว่าพลาดตรงขั้นไหน แทนที่จะต้องรื้อทำใหม่ทั้งห้อง

สรุป

ทั้งงานสรุปได้สั้น ๆ ว่า อย่าเชื่อรูปที่เห็น ให้เชื่อสเกลที่วัดได้ของห้องนั้น ๆ แยกแกน x กับ y เสมอ นับจำนวนด้วย BFS ไม่ใช่ด้วยตา วงกลมใช้ sqrt(w*h) blob ที่สัดส่วนแปลกให้กลับไปซูม แล้วปิดท้ายด้วยการ warp แบบกลับมาเป็นสเกลจริงเพื่อตรวจ

ส่วนตัวผมว่าขั้นที่คุ้มที่สุดคือขั้นวัด px ต่อเมตรรายห้อง ใช้เวลาห้องละไม่กี่นาที แต่มันตัดปัญหาที่จะตามมาทีหลังออกไปเกือบทั้งหมด เทียบกับการวาดผิดแล้วมารู้ตอนแบบเสร็จ ถือว่าถูกมาก

ใครจะเอาไปใช้กับแบบชุดอื่นก็ใช้โครงเดิมได้เลย เปลี่ยนแค่ค่าสีกับกรอบผนัง ส่วนแบบที่เป็นภาพสแกนกระดาษจริง ๆ (มีเงา มีเอียง) อันนั้นต้องเพิ่มขั้น deskew เข้าไปก่อน ผมยังไม่ได้ทำเป็นชุดเก็บไว้ ตรงไหนผมอธิบายตกหล่นหรือพิมพ์ผิด เดี๋ยวกลับมาแก้เพิ่มให้ภายหลังครับ