מודול 8: ניתוח קווי מתאר — ניתוח פורייה אליפטי

שבוע 8 — קווי מתאר ותקופת הברונזה

כשציוני דרך אינם מספיקים

ניתוח המטבעות שלנו השתמש בציוני דרך בדידים — נקודות ניתנות לזיהוי על הדיוקן. אך לחפצים ארכאולוגיים רבים אין נקודות ציון ברורות. חשבו על גרזן ברונזה עם כנפיים: קו מתארו הוא עקומה חלקה ללא נקודות “הומולוגיות” ביולוגית או ארכאולוגית לאורך הקצה.

עבור חפצים כאלה, אנחנו משתמשים בניתוח קווי מתאר — והשיטה הנפוצה ביותר היא ניתוח פורייה אליפטי (EFA).

הרעיון המרכזי: פירוק צורה

EFA מתייחס לקו מתאר כפונקציה מתמטית ומפרק אותו לסכום של רכיבים פשוטים יותר — הרמוניות אליפטיות — בדיוק כשם שאקורד מוזיקלי מורכב ניתן לפירוק לגלי סינוס בודדים (ניתוח פורייה).

האנלוגיה המוזיקלית

תו מוזיקלי הוא גל מורכב
ניתוח פורייה מפרק אותו לגלי סינוס טהורים בתדרים שונים
תדרים נמוכים = האופי הכולל של הצליל
תדרים גבוהים = פרטים עדינים ומרקם

לקווי מתאר: - הרמוניה 1 = הצורה האליפטית הכללית (מאורך לעומת עגול) - הרמוניות 2–4 = תכונות פרופורציה עיקריות (היכן נקודת הרוחב הגדולה ביותר? כמה אסימטרית הבסיס?) - הרמוניות גבוהות יותר = פרטים עדינים (שיניים בקצה, בליטות קטנות)

כיצד EFA עובד

קו מתאר מיוצג כרצף סגור של קואורדינטות $(x_i, y_i)$. EFA מפרק את הקואורדינטות $x$ ו-$y$ בנפרד כסדרות פורייה, ואז משלב אותם מחדש כאליפסות.

כל הרמוניה $n$ מייצרת 4 מקדמים: $(a_n, b_n, c_n, d_n)$.

עבור $N$ הרמוניות, מקבלים $4N$ מקדמים בסך הכל. מקדמים אלה הם הצורה — ניתן להשתמש בהם ישירות ב-PCA בדיוק כמו קואורדינטות ציוני דרך.

נרמול

מקדמי EFA חייבים להיות מנורמלים להסרת אפקטי גודל, מיקום, סיבוב ונקודת התחלה. לאחר נרמול, נותר רק מידע צורה.

בחירת N (מספר הרמוניות)

יותר הרמוניות = יותר פרטים נלכדים, אבל גם יותר רעש. הקריטריון הסטנדרטי:

בחרו את ה-$N$ המינימלי שמסביר ≥ 99% מהשונות בצורת קו המתאר במדגם שלכם.

עבור גרזנים כשלנו, 8–12 הרמוניות בדרך כלל מספיקות.

שלב אחר שלב: EFA ב-Colab

פתחו את מחברת שבוע 8

→ מחברת Colab — שבוע 8

שלב 1: טעינת ועיבוד מקדם של תמונת גרזן אחת

!pip install pyefd scikit-image python-bidi -q

import numpy as np
import matplotlib.pyplot as plt
from bidi.algorithm import get_display
rtl = get_display  # תיקון כיוון טקסט עברי בגרפים של matplotlib
from skimage import io, color, filters, measure
import pyefd

def extract_outline(image_path, threshold_method='otsu'):
    """
    טוען תמונת גרזן וחולץ את קו המתאר הגדול ביותר.
    מחזיר: מערך של קואורדינטות (x, y)
    """
    img = io.imread(image_path, as_gray=True)

    # סף לבינארי
    thresh = filters.threshold_otsu(img)
    binary = img > thresh  # רקע = True, חפץ = False

    # מציאת קווי מתאר של החפץ
    contours = measure.find_contours(~binary, 0.5)

    # שמירת הקו הגדול ביותר (= הגרזן)
    largest = max(contours, key=len)
    return largest  # צורה: (n_נקודות, 2)

# בדיקה על גרזן אחד
outline = extract_outline("data/flanged_axes/G3/ALT.tif")
print(f"לקו המתאר {len(outline)} נקודות")

# ציור
plt.figure(figsize=(6, 8))
plt.plot(outline[:, 1], -outline[:, 0], 'b-', linewidth=1)
plt.fill(outline[:, 1], -outline[:, 0], alpha=0.2)
plt.title(rtl("קו מתאר גרזן ALT"))
plt.axis('equal')
plt.show()

שלב 2: התאמת EFA

# התאמת EFA עם 12 הרמוניות, מנורמל
coeffs = pyefd.elliptic_fourier_descriptors(
    outline, order=12, normalize=True
)
print(f"צורת מקדמי EFA: {coeffs.shape}")  # (12, 4)
# כל שורה = הרמוניה אחת: (a_n, b_n, c_n, d_n)

שלב 3: שחזור הצורה מ-EFA (בדיקת שפיות)

# שחזור קו המתאר ממקדמי EFA (אמור להתאים למקור)
locus = pyefd.reconstruct_contour(coeffs, num_points=300)

plt.figure(figsize=(6, 8))
plt.plot(outline[:, 1], -outline[:, 0], 'gray', linewidth=2, alpha=0.5, label=rtl('מקורי'))
plt.plot(locus[:, 0], -locus[:, 1], 'r-', linewidth=2, label=rtl('שחזור EFA'))
plt.legend()
plt.axis('equal')
plt.title(rtl("שחזור EFA (12 הרמוניות)"))
plt.show()

שלב 4: עיבוד כל הגרזנים ויצירת מטריצה

import os

def process_all_axes(folder_g3, folder_g4):
    """טוען את כל תמונות הגרזן, חולץ קווי מתאר ומתאים EFA."""
    all_coeffs = []
    all_labels = []
    all_names = []

    for folder, group in [(folder_g3, 'G3'), (folder_g4, 'G4')]:
        for fname in sorted(os.listdir(folder)):
            if fname.endswith('.tif'):
                path = os.path.join(folder, fname)
                try:
                    outline = extract_outline(path)
                    coeffs = pyefd.elliptic_fourier_descriptors(
                        outline, order=12, normalize=True
                    )
                    # שיטוח: 12 הרמוניות × 4 מקדמים = 48 ערכים
                    all_coeffs.append(coeffs.flatten())
                    all_labels.append(group)
                    all_names.append(fname.replace('.tif', ''))
                except Exception as e:
                    print(f"דולג {fname}: {e}")

    return np.array(all_coeffs), np.array(all_labels), np.array(all_names)

coeff_matrix, labels, names = process_all_axes(
    "data/flanged_axes/G3/",
    "data/flanged_axes/G4/"
)
print(f"צורת המטריצה: {coeff_matrix.shape}")
# צפוי: (n_גרזנים, 48)
print(f"קבוצות: {dict(zip(*np.unique(labels, return_counts=True)))}")

שלב 5: PCA על מקדמי EFA

from sklearn.decomposition import PCA
import pandas as pd

pca = PCA()
scores = pca.fit_transform(coeff_matrix)
var_exp = pca.explained_variance_ratio_ * 100

df = pd.DataFrame({
    'PC1': scores[:, 0],
    'PC2': scores[:, 1],
    'קבוצה': labels,
    'שם': names
})

# ציור
fig, ax = plt.subplots(figsize=(9, 7))
colors = {'G3': '#1976D2', 'G4': '#C62828'}

for group, grp_df in df.groupby('קבוצה'):
    ax.scatter(grp_df['PC1'], grp_df['PC2'],
               color=colors[group], label=group,
               s=60, alpha=0.85, edgecolor='white')

ax.set_xlabel(f"PC1 ({var_exp[0]:.1f}%)")
ax.set_ylabel(f"PC2 ({var_exp[1]:.1f}%)")
ax.set_title(rtl("מרחב הצורות EFA: גרזני ברונזה עם כנפיים"))
ax.legend(title=rtl("קבוצה טיפולוגית"))
plt.tight_layout()
plt.show()

קריאת גרף PCA של EFA

בניגוד ל-PCA של ציוני דרך, צירי PCA של EFA אינם ממפים ישירות לתזוזות בציוני דרך ספציפיים. במקום זאת, אנחנו מדמיינים את קווי המתאר המשוחזרים בקצוות של כל רכיב עיקרי:

# שחזור צורה ממוצעת וצורות ±2SD לאורך PC1
mean_coeffs = pca.mean_.reshape(12, 4)
pc1_vector  = pca.components_[0].reshape(12, 4)
std1        = np.sqrt(pca.explained_variance_[0])

fig, axes = plt.subplots(1, 3, figsize=(15, 7))
for ax, scale, title in zip(axes, [-2, 0, 2], ['PC1 מינימום (−2SD)', 'צורה ממוצעת', 'PC1 מקסימום (+2SD)']):
    shape_coeffs = mean_coeffs + scale * std1 * pc1_vector
    locus = pyefd.reconstruct_contour(shape_coeffs, num_points=400)
    ax.fill(locus[:, 0], locus[:, 1], alpha=0.3, color='steelblue')
    ax.plot(locus[:, 0], locus[:, 1], 'steelblue', linewidth=2)
    ax.set_title(title)
    ax.set_aspect('equal')
    ax.axis('off')

plt.suptitle("שונות צורה לאורך PC1 — גרזני ברונזה", fontsize=13)
plt.tight_layout()
plt.show()

מה לחפש: האם הצורות המשוחזרות בקצות PC1 מציגות את אותן תכונות מורפולוגיות שמבדילות בין G3 ל-G4 בטיפולוגיה המסורתית?

מטלה 4 מחולקת השבוע

מטלה 4: EFA לגרזני ברונזה — מועד הגשה: שבוע 9

EFA + PCA מלא של גרזנים G3 לעומת G4, עם בדיקה סטטיסטית ומפה גיאוגרפית.

→ מטלה 4

→ מודול 9: מקרה בוחן II — גרזני ברונזה ורשתות אומנות

--- title: "מודול 8: ניתוח קווי מתאר — ניתוח פורייה אליפטי" subtitle: "שבוע 8 — קווי מתאר ותקופת הברונזה" lang: he --- ## כשציוני דרך אינם מספיקים ניתוח המטבעות שלנו השתמש בציוני דרך בדידים — נקודות ניתנות לזיהוי על הדיוקן. אך לחפצים ארכאולוגיים רבים אין נקודות ציון ברורות. חשבו על גרזן ברונזה עם כנפיים: קו מתארו הוא עקומה חלקה ללא נקודות "הומולוגיות" ביולוגית או ארכאולוגית לאורך הקצה. עבור חפצים כאלה, אנחנו משתמשים ב**ניתוח קווי מתאר** — והשיטה הנפוצה ביותר היא **ניתוח פורייה אליפטי (EFA)**. --- ## הרעיון המרכזי: פירוק צורה EFA מתייחס לקו מתאר כפונקציה מתמטית ומפרק אותו לסכום של רכיבים פשוטים יותר — **הרמוניות אליפטיות** — בדיוק כשם שאקורד מוזיקלי מורכב ניתן לפירוק לגלי סינוס בודדים (ניתוח פורייה). ### האנלוגיה המוזיקלית - תו מוזיקלי הוא גל מורכב - ניתוח פורייה מפרק אותו לגלי סינוס טהורים בתדרים שונים - תדרים נמוכים = האופי הכולל של הצליל - תדרים גבוהים = פרטים עדינים ומרקם לקווי מתאר: - **הרמוניה 1** = הצורה האליפטית הכללית (מאורך לעומת עגול) - **הרמוניות 2–4** = תכונות פרופורציה עיקריות (היכן נקודת הרוחב הגדולה ביותר? כמה אסימטרית הבסיס?) - **הרמוניות גבוהות יותר** = פרטים עדינים (שיניים בקצה, בליטות קטנות) --- ## כיצד EFA עובד קו מתאר מיוצג כרצף סגור של קואורדינטות $(x_i, y_i)$. EFA מפרק את הקואורדינטות $x$ ו-$y$ בנפרד כסדרות פורייה, ואז משלב אותם מחדש כאליפסות. כל הרמוניה $n$ מייצרת **4 מקדמים**: $(a_n, b_n, c_n, d_n)$. עבור $N$ הרמוניות, מקבלים $4N$ מקדמים בסך הכל. מקדמים אלה *הם* הצורה — ניתן להשתמש בהם ישירות ב-PCA בדיוק כמו קואורדינטות ציוני דרך. ### נרמול מקדמי EFA חייבים להיות מנורמלים להסרת אפקטי גודל, מיקום, סיבוב ונקודת התחלה. לאחר נרמול, נותר רק מידע צורה. ### בחירת N (מספר הרמוניות) יותר הרמוניות = יותר פרטים נלכדים, אבל גם יותר רעש. הקריטריון הסטנדרטי: > בחרו את ה-$N$ המינימלי שמסביר ≥ 99% מהשונות בצורת קו המתאר במדגם שלכם. עבור גרזנים כשלנו, 8–12 הרמוניות בדרך כלל מספיקות. --- ## שלב אחר שלב: EFA ב-Colab ::: {.callout-important} ## פתחו את מחברת שבוע 8 [→ מחברת Colab — שבוע 8](https://colab.research.google.com/github/shaigordin/comparch/blob/2026/morphometrics/notebooks/week-08.ipynb) ::: ### שלב 1: טעינת ועיבוד מקדם של תמונת גרזן אחת ```python !pip install pyefd scikit-image python-bidi -q import numpy as np import matplotlib.pyplot as plt from bidi.algorithm import get_display rtl = get_display # תיקון כיוון טקסט עברי בגרפים של matplotlib from skimage import io, color, filters, measure import pyefd def extract_outline(image_path, threshold_method='otsu'): """ טוען תמונת גרזן וחולץ את קו המתאר הגדול ביותר. מחזיר: מערך של קואורדינטות (x, y) """ img = io.imread(image_path, as_gray=True) # סף לבינארי thresh = filters.threshold_otsu(img) binary = img > thresh # רקע = True, חפץ = False # מציאת קווי מתאר של החפץ contours = measure.find_contours(~binary, 0.5) # שמירת הקו הגדול ביותר (= הגרזן) largest = max(contours, key=len) return largest # צורה: (n_נקודות, 2) # בדיקה על גרזן אחד outline = extract_outline("data/flanged_axes/G3/ALT.tif") print(f"לקו המתאר {len(outline)} נקודות") # ציור plt.figure(figsize=(6, 8)) plt.plot(outline[:, 1], -outline[:, 0], 'b-', linewidth=1) plt.fill(outline[:, 1], -outline[:, 0], alpha=0.2) plt.title(rtl("קו מתאר גרזן ALT")) plt.axis('equal') plt.show() ``` ### שלב 2: התאמת EFA ```python # התאמת EFA עם 12 הרמוניות, מנורמל coeffs = pyefd.elliptic_fourier_descriptors( outline, order=12, normalize=True ) print(f"צורת מקדמי EFA: {coeffs.shape}") # (12, 4) # כל שורה = הרמוניה אחת: (a_n, b_n, c_n, d_n) ``` ### שלב 3: שחזור הצורה מ-EFA (בדיקת שפיות) ```python # שחזור קו המתאר ממקדמי EFA (אמור להתאים למקור) locus = pyefd.reconstruct_contour(coeffs, num_points=300) plt.figure(figsize=(6, 8)) plt.plot(outline[:, 1], -outline[:, 0], 'gray', linewidth=2, alpha=0.5, label=rtl('מקורי')) plt.plot(locus[:, 0], -locus[:, 1], 'r-', linewidth=2, label=rtl('שחזור EFA')) plt.legend() plt.axis('equal') plt.title(rtl("שחזור EFA (12 הרמוניות)")) plt.show() ``` ### שלב 4: עיבוד כל הגרזנים ויצירת מטריצה ```python import os def process_all_axes(folder_g3, folder_g4): """טוען את כל תמונות הגרזן, חולץ קווי מתאר ומתאים EFA.""" all_coeffs = [] all_labels = [] all_names = [] for folder, group in [(folder_g3, 'G3'), (folder_g4, 'G4')]: for fname in sorted(os.listdir(folder)): if fname.endswith('.tif'): path = os.path.join(folder, fname) try: outline = extract_outline(path) coeffs = pyefd.elliptic_fourier_descriptors( outline, order=12, normalize=True ) # שיטוח: 12 הרמוניות × 4 מקדמים = 48 ערכים all_coeffs.append(coeffs.flatten()) all_labels.append(group) all_names.append(fname.replace('.tif', '')) except Exception as e: print(f"דולג {fname}: {e}") return np.array(all_coeffs), np.array(all_labels), np.array(all_names) coeff_matrix, labels, names = process_all_axes( "data/flanged_axes/G3/", "data/flanged_axes/G4/" ) print(f"צורת המטריצה: {coeff_matrix.shape}") # צפוי: (n_גרזנים, 48) print(f"קבוצות: {dict(zip(*np.unique(labels, return_counts=True)))}") ``` ### שלב 5: PCA על מקדמי EFA ```python from sklearn.decomposition import PCA import pandas as pd pca = PCA() scores = pca.fit_transform(coeff_matrix) var_exp = pca.explained_variance_ratio_ * 100 df = pd.DataFrame({ 'PC1': scores[:, 0], 'PC2': scores[:, 1], 'קבוצה': labels, 'שם': names }) # ציור fig, ax = plt.subplots(figsize=(9, 7)) colors = {'G3': '#1976D2', 'G4': '#C62828'} for group, grp_df in df.groupby('קבוצה'): ax.scatter(grp_df['PC1'], grp_df['PC2'], color=colors[group], label=group, s=60, alpha=0.85, edgecolor='white') ax.set_xlabel(f"PC1 ({var_exp[0]:.1f}%)") ax.set_ylabel(f"PC2 ({var_exp[1]:.1f}%)") ax.set_title(rtl("מרחב הצורות EFA: גרזני ברונזה עם כנפיים")) ax.legend(title=rtl("קבוצה טיפולוגית")) plt.tight_layout() plt.show() ``` --- ## קריאת גרף PCA של EFA בניגוד ל-PCA של ציוני דרך, צירי PCA של EFA אינם ממפים ישירות לתזוזות בציוני דרך ספציפיים. במקום זאת, אנחנו מדמיינים את *קווי המתאר המשוחזרים* בקצוות של כל רכיב עיקרי: ```python # שחזור צורה ממוצעת וצורות ±2SD לאורך PC1 mean_coeffs = pca.mean_.reshape(12, 4) pc1_vector = pca.components_[0].reshape(12, 4) std1 = np.sqrt(pca.explained_variance_[0]) fig, axes = plt.subplots(1, 3, figsize=(15, 7)) for ax, scale, title in zip(axes, [-2, 0, 2], ['PC1 מינימום (−2SD)', 'צורה ממוצעת', 'PC1 מקסימום (+2SD)']): shape_coeffs = mean_coeffs + scale * std1 * pc1_vector locus = pyefd.reconstruct_contour(shape_coeffs, num_points=400) ax.fill(locus[:, 0], locus[:, 1], alpha=0.3, color='steelblue') ax.plot(locus[:, 0], locus[:, 1], 'steelblue', linewidth=2) ax.set_title(title) ax.set_aspect('equal') ax.axis('off') plt.suptitle("שונות צורה לאורך PC1 — גרזני ברונזה", fontsize=13) plt.tight_layout() plt.show() ``` **מה לחפש:** האם הצורות המשוחזרות בקצות PC1 מציגות את אותן תכונות מורפולוגיות שמבדילות בין G3 ל-G4 בטיפולוגיה המסורתית? --- ## מטלה 4 מחולקת השבוע ::: {.callout-note} ## מטלה 4: EFA לגרזני ברונזה — מועד הגשה: שבוע 9 EFA + PCA מלא של גרזנים G3 לעומת G4, עם בדיקה סטטיסטית ומפה גיאוגרפית. [→ מטלה 4](../assignments/04-axes.qmd) ::: [→ מודול 9: מקרה בוחן II — גרזני ברונזה ורשתות אומנות](09-statistics.qmd)