پرش به محتویات

یافتن نزدیک‌ترین جفت نقاط

صورت مسئله

$n$ نقطه در صفحه داده شده است. هر نقطه $p_i$ با مختصات $(x_i,y_i)$ خود تعریف می‌شود. هدف، یافتن دو نقطه از میان آن‌هاست به طوری که فاصله بینشان کمینه باشد:

$$ \min_{\scriptstyle i, j=0 \ldots n-1,\atop \scriptstyle i \neq j } \rho (p_i, p_j). $$

ما از فاصله اقلیدسی معمول استفاده می‌کنیم:

$$ \rho (p_i,p_j) = \sqrt{(x_i-x_j)^2 + (y_i-y_j)^2} .$$

الگوریتم بدیهی - یعنی پیمایش تمام جفت‌ها و محاسبه فاصله برای هر کدام - در زمان $O(n^2)$ اجرا می‌شود.

در ادامه، الگوریتمی با زمان اجرای $O(n \log n)$ شرح داده می‌شود. این الگوریتم توسط Shamos و Hoey در سال ۱۹۷۵ ارائه شد. (منبع: یادداشت‌های فصل ۵ کتاب Algorithm Design نوشته Kleinberg & Tardos، همچنین اینجا را ببینید). Preparata و Shamos همچنین نشان دادند که این الگوریتم در مدل درخت تصمیم (decision tree model) بهینه است.

الگوریتم

ما الگوریتمی را بر اساس طرح کلی الگوریتم‌های تقسیم و حل (divide-and-conquer) می‌سازیم: الگوریتم به صورت یک تابع بازگشتی طراحی می‌شود که مجموعه‌ای از نقاط را به عنوان ورودی می‌گیرد؛ این تابع بازگشتی مجموعه را به دو نیمه تقسیم کرده، خود را به صورت بازگشتی روی هر نیمه فراخوانی می‌کند و سپس عملیاتی را برای ترکیب جواب‌ها انجام می‌دهد. عملیات ترکیب شامل شناسایی مواردی است که یک نقطه از جفت بهینه در یک نیمه و نقطه دیگر در نیمه دیگر قرار گرفته باشد (در این حالت، فراخوانی‌های بازگشتی از هر یک از نیمه‌ها به تنهایی نمی‌توانند این جفت را تشخیص دهند). مشکل اصلی، مانند همیشه در الگوریتم‌های تقسیم و حل، در پیاده‌سازی مؤثر مرحله ادغام نهفته است. اگر مجموعه‌ای از $n$ نقطه به تابع بازگشتی داده شود، آنگاه مرحله ادغام نباید بیشتر از $O(n)$ طول بکشد، در این صورت پیچیدگی کل الگوریتم $T(n)$ از معادله زیر به دست می‌آید:

$$T(n) = 2T(n/2) + O(n).$$

همانطور که می‌دانیم، حل این معادله $T(n) = O(n \log n)$ است.

بنابراین، به ساخت الگوریتم می‌پردازیم. برای رسیدن به یک پیاده‌سازی کارآمد برای مرحله ادغام در آینده، مجموعه نقاط را بر اساس مختصات $x$ آن‌ها به دو زیرمجموعه تقسیم می‌کنیم: در واقع، یک خط عمودی رسم می‌کنیم که مجموعه نقاط را به دو زیرمجموعه با اندازه‌های تقریباً برابر تقسیم کند. انجام چنین تقسیمی به این صورت راحت است: نقاط را به صورت استاندارد به عنوان جفت اعداد مرتب می‌کنیم، یعنی:

$$p_i < p_j \Longleftrightarrow (x_i < x_j) \lor \Big(\left(x_i = x_j\right) \wedge \left(y_i < y_j \right) \Big) $$

سپس نقطه میانی پس از مرتب‌سازی یعنی $p_m (m = \lfloor n/2 \rfloor)$ را در نظر می‌گیریم و تمام نقاط قبل از آن به همراه خود $p_m$ را به نیمه اول، و تمام نقاط بعد از آن را به نیمه دوم اختصاص می‌دهیم:

$$A_1 = \{p_i \ | \ i = 0 \ldots m \}$$
$$A_2 = \{p_i \ | \ i = m + 1 \ldots n-1 \}.$$

حالا با فراخوانی بازگشتی روی هر یک از مجموعه‌های $A_1$ و $A_2$، جواب‌های $h_1$ و $h_2$ را برای هر نیمه پیدا می‌کنیم و بهترین آن‌ها را انتخاب می‌کنیم: $h = \min(h_1, h_2)$.

اکنون باید مرحله ادغام را انجام دهیم، یعنی سعی می‌کنیم جفت نقاطی را پیدا کنیم که فاصله بین آن‌ها کمتر از $h$ باشد و یک نقطه در $A_1$ و دیگری در $A_2$ قرار داشته باشد. بدیهی است که کافی است تنها نقاطی را در نظر بگیریم که از خط عمودی فاصله‌ای کمتر از $h$ دارند، یعنی مجموعه $B$ از نقاطی که در این مرحله بررسی می‌شوند برابر است با:

$$B = \{ p_i\ | \ | x_i - x_m\ | < h \}.$$

برای هر نقطه در مجموعه $B$، سعی می‌کنیم نقاطی را پیدا کنیم که از آن به $h$ نزدیک‌تر باشند. به عنوان مثال، کافی است تنها نقاطی را در نظر بگیریم که مختصات $y$ آن‌ها حداکثر $h$ تفاوت داشته باشد. علاوه بر این، در نظر گرفتن نقاطی که مختصات $y$ آن‌ها بزرگتر از مختصات $y$ نقطه فعلی است، بی‌معنی است. بنابراین، برای هر نقطه $p_i$ مجموعه نقاط مورد بررسی $C(p_i)$ را به صورت زیر تعریف می‌کنیم:

$$C(p_i) = \{ p_j\ |\ p_j \in B,\ \ y_i - h < y_j \le y_i \}.$$

اگر نقاط مجموعه $B$ را بر اساس مختصات $y$ مرتب کنیم، پیدا کردن $C(p_i)$ بسیار آسان خواهد بود: این‌ها چند نقطه متوالی قبل از نقطه $p_i$ هستند.

بنابراین، با نمادگذاری جدید، مرحله ادغام به این صورت است: مجموعه $B$ را بسازید، نقاط موجود در آن را بر اساس مختصات $y$ مرتب کنید، سپس برای هر نقطه $p_i \in B$ تمام نقاط $p_j \in C(p_i)$ را در نظر بگیرید و برای هر جفت $(p_i,p_j)$ فاصله را محاسبه کرده و با بهترین فاصله فعلی مقایسه کنید.

در نگاه اول، این هنوز یک الگوریتم غیربهینه به نظر می‌رسد: به نظر می‌رسد که اندازه مجموعه‌های $C(p_i)$ از مرتبه $n$ خواهد بود و پیچیدگی مورد نظر به دست نخواهد آمد. با این حال، به طرز شگفت‌آوری، می‌توان ثابت کرد که اندازه هر یک از مجموعه‌های $C(p_i)$ یک مقدار $O(1)$ است، یعنی صرف‌نظر از خود نقاط، از یک ثابت کوچک تجاوز نمی‌کند. اثبات این واقعیت در بخش بعدی آورده شده است.

در نهایت، به مرتب‌سازی‌هایی که الگوریتم فوق شامل می‌شود توجه می‌کنیم: اول، مرتب‌سازی بر اساس جفت‌های $(x, y)$ و دوم، مرتب‌سازی عناصر مجموعه $B$ بر اساس $y$. در واقع، هر دوی این مرتب‌سازی‌ها را می‌توان از داخل تابع بازگشتی حذف کرد (در غیر این صورت به تخمین $O(n)$ برای مرحله ادغام نمی‌رسیدیم و پیچیدگی کلی الگوریتم $O(n \log^2 n)$ می‌شد). خلاص شدن از مرتب‌سازی اول آسان است — کافی است این مرتب‌سازی را قبل از شروع بازگشت انجام دهیم: چرا که خود عناصر در داخل بازگشت تغییر نمی‌کنند، بنابراین نیازی به مرتب‌سازی مجدد نیست. انجام مرتب‌سازی دوم کمی دشوارتر است و انجام آن از قبل کارساز نخواهد بود. اما، با به یاد آوردن مرتب‌سازی ادغامی (merge sort) که آن هم بر اساس اصل تقسیم و حل کار می‌کند، می‌توانیم این مرتب‌سازی را به سادگی در بازگشت خود جای دهیم. اجازه دهید تابع بازگشتی، با گرفتن یک مجموعه از نقاط (که همانطور که به یاد داریم، بر اساس جفت‌های $(x, y)$ مرتب شده‌اند)، همان مجموعه را برگرداند، اما این بار مرتب شده بر اساس مختصات $y$. برای انجام این کار، کافی است دو نتیجه بازگشتی را با هم ادغام کنیم (در زمان $O(n)$). این کار منجر به مجموعه‌ای مرتب شده بر اساس مختصات $y$ خواهد شد.

ارزیابی پیچیدگی زمانی

برای نشان دادن اینکه الگوریتم فوق در واقع در زمان $O(n \log n)$ اجرا می‌شود، باید واقعیت زیر را ثابت کنیم: $|C(p_i)| = O(1)$.

بنابراین، یک نقطه $p_i$ را در نظر می‌گیریم؛ به یاد بیاورید که مجموعه $C(p_i)$ مجموعه‌ای از نقاط است که مختصات $y$ آن‌ها در بازه $[y_i-h; y_i]$ قرار دارد و علاوه بر این، در امتداد مختصات $x$، خود نقطه $p_i$ و تمام نقاط مجموعه $C(p_i)$ در نواری به عرض $2h$ قرار دارند. به عبارت دیگر، نقاطی که ما در نظر می‌گیریم، یعنی $p_i$ و $C(p_i)$، در یک مستطیل به ابعاد $2h \times h$ قرار دارند.

وظیفه ما تخمین حداکثر تعداد نقاطی است که می‌توانند در این مستطیل $2h \times h$ قرار بگیرند؛ بنابراین، حداکثر اندازه مجموعه $C(p_i)$ را تخمین می‌زنیم. در عین حال، هنگام ارزیابی، نباید فراموش کنیم که ممکن است نقاط تکراری وجود داشته باشند.

به یاد داشته باشید که $h$ از نتایج دو فراخوانی بازگشتی - روی مجموعه‌های $A_1$ و $A_2$ - به دست آمده است، و $A_1$ شامل نقاط سمت چپ خط تقسیم و بخشی از نقاط روی آن است، و $A_2$ شامل نقاط باقی‌مانده روی خط تقسیم و نقاط سمت راست آن است. برای هر جفت نقطه از $A_1$ و همچنین از $A_2$، فاصله نمی‌تواند کمتر از $h$ باشد — در غیر این صورت به معنای عملکرد نادرست تابع بازگشتی خواهد بود.

برای تخمین حداکثر تعداد نقاط در مستطیل $2h \times h$، آن را به دو مربع $h \times h$ تقسیم می‌کنیم. مربع اول شامل تمام نقاط $C(p_i) \cap A_1$ و مربع دوم شامل بقیه، یعنی $C(p_i) \cap A_2$ است. از ملاحظات فوق نتیجه می‌شود که در هر یک از این مربع‌ها فاصله بین هر دو نقطه حداقل $h$ است.

نشان می‌دهیم که در هر مربع حداکثر چهار نقطه وجود دارد. به عنوان مثال، این کار را می‌توان به صورت زیر انجام داد: مربع را به ۴ زیرمربع با اضلاع $h/2$ تقسیم کنید. آنگاه در هر یک از این زیرمربع‌ها نمی‌تواند بیش از یک نقطه وجود داشته باشد (زیرا حتی قطر آن برابر با $h / \sqrt{2}$ است که کمتر از $h$ است). بنابراین، در کل مربع نمی‌تواند بیش از ۴ نقطه وجود داشته باشد.

بنابراین، ما ثابت کردیم که در یک مستطیل $2h \times h$ نمی‌تواند بیش از $4 \cdot 2 = 8$ نقطه وجود داشته باشد، و در نتیجه، اندازه مجموعه $C(p_i)$ نمی‌تواند از ۷ تجاوز کند، که همان چیزی بود که می‌خواستیم.

پیاده‌سازی

یک ساختار داده برای ذخیره یک نقطه (مختصات و شماره آن) و عملگرهای مقایسه مورد نیاز برای دو نوع مرتب‌سازی معرفی می‌کنیم:

struct pt {
    int x, y, id;
};

struct cmp_x {
    bool operator()(const pt & a, const pt & b) const {
        return a.x < b.x || (a.x == b.x && a.y < b.y);
    }
};

struct cmp_y {
    bool operator()(const pt & a, const pt & b) const {
        return a.y < b.y;
    }
};

int n;
vector<pt> a;

برای پیاده‌سازی راحت‌تر بازگشت، یک تابع کمکی upd_ans() معرفی می‌کنیم که فاصله بین دو نقطه را محاسبه کرده و بررسی می‌کند که آیا از جواب فعلی بهتر است یا خیر:

double mindist;
pair<int, int> best_pair;

void upd_ans(const pt & a, const pt & b) {
    double dist = sqrt((a.x - b.x)*(a.x - b.x) + (a.y - b.y)*(a.y - b.y));
    if (dist < mindist) {
        mindist = dist;
        best_pair = {a.id, b.id};
    }
}

در نهایت، پیاده‌سازی خود بازگشت. فرض بر این است که قبل از فراخوانی آن، آرایه a[] از قبل بر اساس مختصات $x$ مرتب شده است. در بازگشت فقط دو اشاره‌گر l و r را پاس می‌دهیم که نشان می‌دهد باید به دنبال پاسخ برای a[l ... r) بگردد. اگر فاصله بین r و l خیلی کم باشد، بازگشت باید متوقف شود و یک الگوریتم بدیهی برای یافتن نزدیک‌ترین جفت اجرا شود و سپس زیرآرایه بر اساس مختصات y مرتب شود.

برای ادغام دو مجموعه از نقاط دریافت شده از فراخوانی‌های بازگشتی به یک مجموعه واحد (مرتب شده بر اساس مختصات y)، از تابع استاندارد merge() از STL استفاده می‌کنیم و یک بافر کمکی t[] (یکی برای تمام فراخوانی‌های بازگشتی) ایجاد می‌کنیم. (استفاده از inplace_merge() عملی نیست زیرا به طور کلی در زمان خطی کار نمی‌کند.)

در نهایت، مجموعه $B$ در همان آرایه t ذخیره می‌شود.

vector<pt> t;

void rec(int l, int r) {
    if (r - l <= 3) {
        for (int i = l; i < r; ++i) {
            for (int j = i + 1; j < r; ++j) {
                upd_ans(a[i], a[j]);
            }
        }
        sort(a.begin() + l, a.begin() + r, cmp_y());
        return;
    }

    int m = (l + r) >> 1;
    int midx = a[m].x;
    rec(l, m);
    rec(m, r);

    merge(a.begin() + l, a.begin() + m, a.begin() + m, a.begin() + r, t.begin(), cmp_y());
    copy(t.begin(), t.begin() + r - l, a.begin() + l);

    int tsz = 0;
    for (int i = l; i < r; ++i) {
        if (abs(a[i].x - midx) < mindist) {
            for (int j = tsz - 1; j >= 0 && a[i].y - t[j].y < mindist; --j)
                upd_ans(a[i], t[j]);
            t[tsz++] = a[i];
        }
    }
}

ضمناً، اگر تمام مختصات صحیح باشند، می‌توانید در طول بازگشت به مقادیر کسری نروید و در mindist مربع کمترین فاصله را ذخیره کنید.

در برنامه اصلی، تابع بازگشتی باید به این صورت فراخوانی شود:

t.resize(n);
sort(a.begin(), a.end(), cmp_x());
mindist = 1E20;
rec(0, n);

تعمیم: یافتن مثلثی با کمترین محیط

الگوریتم توصیف شده در بالا به طرز جالبی به این مسئله تعمیم می‌یابد: از میان مجموعه‌ای از نقاط داده شده، سه نقطه مختلف را به گونه‌ای انتخاب کنید که مجموع فواصل دو به دوی بین آن‌ها کمترین مقدار ممکن باشد.

در واقع، برای حل این مسئله، الگوریتم به همان صورت باقی می‌ماند: صفحه را با یک خط عمودی به دو نیمه تقسیم می‌کنیم، راه حل را به صورت بازگشتی روی هر دو نیمه فراخوانی می‌کنیم، کمترین محیط minper را از میان محیط‌های یافت شده انتخاب می‌کنیم، نواری به ضخامت minper / 2 می‌سازیم و تمام مثلث‌هایی که می‌توانند پاسخ را بهبود بخشند، بررسی می‌کنیم. (توجه داشته باشید که مثلثی با محیط $\le minper$ دارای بزرگترین ضلع $\le minper / 2$ است.)

مسائل تمرینی