From octave-sources-request at bevo dot che dot wisc dot edu Tue Dec 30 22:20:01 2003 Subject: miniball - Smallest enclosing ball From: Harald Anlauf To: octave-sources at bevo dot che dot wisc dot edu Date: Sun, 28 Dec 2003 11:31:30 -0600 This is a MIME encoded message. --NHgP2=-bQA=_U4Sw=?9SSKsrLnwKdFMoWqGDjv Content-Type: text/plain; charset=iso-8859-1 Hi all, the appended algorithm solves a common problem in computational geometry. It determines the (euclidian) center and radius of the smallest enclosing ball of a set of n points in d dimensions. Example: 50 random points in 5 dimensions octave:1> a=randn(5,50); [c,r]=miniball(a) c = -0.491455 0.475413 -0.685300 0.096564 -0.164695 r = 3.3357 I hereby release the code under the GPL for inclusion in a future release of octave-(forge). Comments are welcome. Have fun, -Harald --NHgP2=-bQA=_U4Sw=?9SSKsrLnwKdFMoWqGDjv Content-Type: application/octet-stream; name="miniball.m" Content-Transfer-Encoding: base64 Content-Disposition: attachment; filename="miniball.m" IyMgQ29weXJpZ2h0IChDKSAyMDAzIEhhcmFsZCBBbmxhdWYKIyMKIyMgVGhpcyBwcm9ncmFt IGlzIGZyZWUgc29mdHdhcmU7IHlvdSBjYW4gcmVkaXN0cmlidXRlIGl0IGFuZC9vciBtb2Rp ZnkKIyMgaXQgdW5kZXIgdGhlIHRlcm1zIG9mIHRoZSBHTlUgR2VuZXJhbCBQdWJsaWMgTGlj ZW5zZSBhcyBwdWJsaXNoZWQgYnkKIyMgdGhlIEZyZWUgU29mdHdhcmUgRm91bmRhdGlvbjsg ZWl0aGVyIHZlcnNpb24gMiBvZiB0aGUgTGljZW5zZSwgb3IKIyMgKGF0IHlvdXIgb3B0aW9u KSBhbnkgbGF0ZXIgdmVyc2lvbi4KIyMKIyMgVGhpcyBwcm9ncmFtIGlzIGRpc3RyaWJ1dGVk IGluIHRoZSBob3BlIHRoYXQgaXQgd2lsbCBiZSB1c2VmdWwsCiMjIGJ1dCBXSVRIT1VUIEFO WSBXQVJSQU5UWTsgd2l0aG91dCBldmVuIHRoZSBpbXBsaWVkIHdhcnJhbnR5IG9mCiMjIE1F UkNIQU5UQUJJTElUWSBvciBGSVRORVNTIEZPUiBBIFBBUlRJQ1VMQVIgUFVSUE9TRS4gIFNl ZSB0aGUKIyMgR05VIEdlbmVyYWwgUHVibGljIExpY2Vuc2UgZm9yIG1vcmUgZGV0YWlscy4K IyMKIyMgWW91IHNob3VsZCBoYXZlIHJlY2VpdmVkIGEgY29weSBvZiB0aGUgR05VIEdlbmVy YWwgUHVibGljIExpY2Vuc2UKIyMgYWxvbmcgd2l0aCB0aGlzIHByb2dyYW07IGlmIG5vdCwg d3JpdGUgdG8gdGhlIEZyZWUgU29mdHdhcmUKIyMgRm91bmRhdGlvbiwgSW5jLiwgNTkgVGVt cGxlIFBsYWNlLCBTdWl0ZSAzMzAsIEJvc3RvbiwgTUEgIDAyMTExLTEzMDcgIFVTQSwKIyMg b3IgZG93bmxvYWQgdGhlIExpY2Vuc2UgdGVybXMgZnJvbSBodHRwOi8vd3d3LmdudS5vcmcv bGljZW5zZXMvZ3BsLnR4dAoKIyMgLSotIHRleGluZm8gLSotCiMjIEBkZWZ0eXBlZm4ge0Z1 bmN0aW9uIEZpbGV9IHtbQHZhcntjZW50ZXJ9LCBAdmFye3JhZGl1c31dID19IG1pbmliYWxs IChAdmFye3BvaW50c30pCiMjIENvbXB1dGUgdGhlIHNtYWxsZXN0IGVuY2xvc2luZyBiYWxs IGZvciBuIHBvaW50cyBvZiBkaW1lbnNpb24gZAojIyBzdG9yZWQgaW4gdGhlIGlucHV0IG1h dHJpeCBAdmFye3BvaW50c30gb2Ygc2l6ZSBbZCwgbl0uICBVcG9uIHJldHVybiwKIyMgQHZh cntjZW50ZXJ9IGFuZCBAdmFye3JhZGl1c30gY29udGFpbiB0aGUgY2VudGVyIGFuZCByYWRp dXMgb2YgdGhlIGJhbGwuCiMjIEBlbmQgZGVmdHlwZWZuCgpmdW5jdGlvbiBbY2VudGVyLHJh ZGl1c10gPSBtaW5pYmFsbCAocG9pbnRzZXQpCgojIyBBdXRob3I6IEhhcmFsZCBBbmxhdWYg PGFubGF1ZiBhdCBoZXAgZG90IHR1LWRhcm1zdGFkdCBkb3QgZGU+CiMjIEZpcnN0IFJlbGVh c2U6IDIwMDMtMTItMjgKIyMgVGhpcyBzb2Z0d2FyZSBpcyBkaXN0cmlidXRlZCB1bmRlciB0 aGUgdGVybXMgb2YgdGhlIEdQTAoKICBbZCxuXSA9IHNpemUocG9pbnRzZXQpOwoKICBpZiAo ZCA9PSAwIHx8IG4gPT0gMCkKICAgIGVycm9yICgnbWluaWJhbGw6IG5vbmNvbmZvcm1pbmcg YXJndW1lbnQ6IHBvaW50c2V0IG11c3QgYmUgYXJyYXkgb2Ygc2l6ZSAoZCxuKScpOwogIGVu ZGlmCgogICUlIEhhbmRsZSB0cml2aWFsIGNhc2VzIHNlcGFyYXRlbHkKICBpZiAobiA9PSAx IHx8IG4gPT0gMikKICAgIFtjZW50ZXIscmFkaXVzX3Nxcl0gPSBtaW5pYmFsbF9iZnNzKHBv aW50c2V0KTsKICAgIHJhZGl1cyA9IHNxcnQocmFkaXVzX3Nxcik7CiAgICByZXR1cm4KICBl bmRpZgoKICAlJSBEZXRlcm1pbmF0aW9uIG9mIHRoZSBzdXBwb3J0IHNldCBmb3IgdGhlIGVu Y2xvc2luZyBiYWxsLCBpbnNwaXJlZAogICUlIGJ5IEfkcnRuZXIncyBbMV0gcGl2b3Rpbmcg dmFyaWFudCBvZiBXZWx6bCdzIFsyXSAibW92ZS10by1mcm9udCIKICAlJSBzY2hlbWUuICBI b3dldmVyLCB0aGUgaW1wbGVtZW50YXRpb24gYmVsb3cgZW1wbG95cyBhIGRpZmZlcmVudAog ICUlIGFsZ29yaXRobSB0byBzZWxlY3QgdGhlIHBvaW50cyBvZiB0aGUgc3VwcG9ydCBzZXQu CiAgJSUKICAlJSBbMV0gQi4gR+RydG5lciwgIkZhc3QgYW5kIFJvYnVzdCBTbWFsbGVzdCBF bmNsb3NpbmcgQmFsbHMiLCBpbjoKICAlJSAgICAgSi4gTmVzZXRyaWwgKGVkLiksIEVTQSc5 OSwgTE5DUyAxNjMsIHBwLiAzMjUtMzM4LCBTcHJpbmdlciwgMTk5OS4KICAlJSBbMl0gRS4g V2VsemwsICJTbWFsbGVzdCBlbmNsb3NpbmcgZGlza3MgKGJhbGxzIGFuZCBlbGxpcHNvaWRz KSIsIGluOgogICUlICAgICBILiBNYXVyZXIgKGVkLiksICJOZXcgUmVzdWx0cyBhbmQgTmV3 IFRyZW5kcyBpbiBDb21wdXRlciBTY2llbmNlIiwKICAlJSAgICAgTGVjdHVyZSBOb3RlcyBp biBDb21wdXRlciBTY2llbmNlIDU1NSAoMTk5MSkgMzU5LTM3MC4KICAlJQogICUlIFRoZSBj dXJyZW50IGNhbmRpZGF0ZXMgZm9yIHRoZSBzdXBwb3J0IHNldCBhcmUgbW92ZWQgdG8KICAl JSBwb2ludHNldCg6LDE6bSksIG0gPD0gZCsxCgogICUlIEluaXRpYWwgY2FuZGlkYXRlIGZv ciBzdXBwb3J0IHNldAogIGN1cl9jZW50ZXIgPSBwb2ludHNldCg6LDEpOwogIGN1cl9yMiA9 IDA7CiAgZV9tYXggPSAwOwogIGlfbWF4ID0gMDsKICBmb3IgaSA9IDI6bgogICAgZSA9IG5v cm0ocG9pbnRzZXQoOixpKS1jdXJfY2VudGVyKV4yIC0gY3VyX3IyOwogICAgaWYgKGUgPiBl X21heCkKICAgICAgZV9tYXggPSBlOwogICAgICBpX21heCA9IGk7CiAgICBlbmRpZgogIGVu ZGZvcgogIGlmIChpX21heCA9PSAwKQogICAgd2FybmluZygibWluaWJhbGw6IGNhbGxlZCB3 aXRoIGRlZ2VuZXJhdGUgc2V0IG9mIHBvaW50cyEiKTsKICAgIGNlbnRlciA9IGN1cl9jZW50 ZXI7CiAgICByYWRpdXMgPSByYWRpdXNfc3FyID0gMDsKICAgIHJldHVybgogIGVuZGlmCgog IGlmIChpX21heCA+PSAwKQogICAgJSUgU3dhcCBwb2ludHMKICAgIG5ld3AgPSBwb2ludHNl dCg6LGlfbWF4KTsKICAgIHBvaW50c2V0KDosaV9tYXgpID0gcG9pbnRzZXQoOiwxKTsKICAg IHBvaW50c2V0KDosMSkgPSBuZXdwOwogIGVuZGlmCiAgJSUgSW5pdGlhbCBzaXplIG9mIHN1 cHBvcnQgc2V0CiAgbSA9IDE7CgogIHRvbCA9IDFlLTE1OwogIGl0ZXIgPSAwOwogICUlIEl0 ZXJhdGlvbiBsaW1pdCBmb3IgZGVnZW5lcmF0ZSBjYXNlcwogIGl0ZXJfbWF4ID0gbWF4KGQr MSxuKTsKICB3aGlsZSAoaXRlciA8IGl0ZXJfbWF4KQogICAgJSUgU2VhcmNoIGZvciBwb2lu dCBvdXRzaWRlIGN1cnJlbnQgYmFsbCB3aXRoICJtYXhpbXVtIGV4Y2VzcyIKICAgIGVfbWF4 ID0gMDsKICAgIGlfbWF4ID0gMDsKICAgIGl0ZXIrKzsKICAgIGZvciBpID0gbSsxOm4KICAg ICAgZSA9IG5vcm0ocG9pbnRzZXQoOixpKS1jdXJfY2VudGVyKV4yIC0gY3VyX3IyOwogICAg ICBpZiAoZSA+IGVfbWF4KQoJZV9tYXggPSBlOwoJaV9tYXggPSBpOwogICAgICBlbmRpZgog ICAgZW5kZm9yCiAgICBpZiAoaV9tYXggPT0gMCB8fCBlX21heCA8IGN1cl9yMip0b2wpCiAg ICAgIGJyZWFrOwogICAgZW5kaWYKICAgICUlIE1vdmUgdG8gZnJvbnQKICAgIGlmIChtID09 IGQrMSkKICAgICAgbmV3bSA9IG07CiAgICBlbHNlCiAgICAgIG5ld20gPSBtKzE7CiAgICBl bmRpZgogICAgbmV3cCA9IHBvaW50c2V0KDosaV9tYXgpOwogICAgaWYgKGlfbWF4IH49IG5l d20pCiAgICAgIHBvaW50c2V0KDosaV9tYXgpID0gcG9pbnRzZXQoOixuZXdtKTsKICAgIGVu ZGlmCiAgICBmb3IgaSA9IG5ld206LTE6MgogICAgICBwb2ludHNldCg6LGkpID0gcG9pbnRz ZXQoOixpLTEpOwogICAgZW5kZm9yCiAgICBwb2ludHNldCg6LDEpID0gbmV3cDsKICAgIG0g PSBuZXdtOwogICAgJSUgRGV0ZXJtaW5lIG5ldyBiYWxsIGNhbmRpZGF0ZQogICAgW2N1cl9j ZW50ZXIsY3VyX3IyLGxhbWJkYV0gPSBtaW5pYmFsbF9iZnNzKHBvaW50c2V0KDosMTptKSk7 CiAgICAlJSBDaGVjayB3aGV0aGVyIHRoZSBiYWxsJ3MgY2VudGVyIGxpZXMgd2l0aGluIHRo ZSBzaW1wbGV4CiAgICB3aGlsZSAoKGxhX21pbiA9IG1pbihsYW1iZGEpKSA8IDApCiAgICAg ICUlIERyb3AgImJhZCIgcG9pbnRzIGZyb20gc3VwcG9ydCBzZXQKICAgICAgbWF4bSA9IG07 CiAgICAgIGZvciBpID0gbTotMToxCglpZiAobGFtYmRhKGkpIDwgMCkKCSAgaWYgKGkgPT0g bWF4bSkKCSAgICBtYXhtLS07CgkgIGVsc2UKCSAgICBuZXdwID0gcG9pbnRzZXQoOixtYXht KTsKCSAgICBwb2ludHNldCg6LG1heG0pID0gcG9pbnRzZXQoOixpKTsKCSAgICBwb2ludHNl dCg6LGkpID0gbmV3cDsKCSAgICBtYXhtLS07CgkgIGVuZGlmCgllbmRpZgogICAgICBlbmRm b3IKICAgICAgbSA9IG1heG07CiAgICAgIFtjdXJfY2VudGVyLGN1cl9yMixsYW1iZGFdID0g bWluaWJhbGxfYmZzcyhwb2ludHNldCg6LDE6bSkpOwogICAgZW5kd2hpbGUKICBlbmR3aGls ZQoKICBjZW50ZXIgPSBjdXJfY2VudGVyOwogIHJhZGl1c19zcXIgPSBjdXJfcjI7CiAgcmFk aXVzID0gc3FydChyYWRpdXNfc3FyKTsKCiAgaWYgKGl0ZXIgPT0gaXRlcl9tYXgpCiAgICB3 YXJuaW5nKCJJdGVyYXRpb24gbGltaXQgcmVhY2hlZCwgcmVzdWx0IG9mIG1pbmliYWxsIG1h eSBiZSBpbmFjY3VyYXRlIik7CiAgZW5kaWYKCmVuZGZ1bmN0aW9uCgojIyBCYWxsIGZyb20g c3VwcG9ydCBzZXQuICBJbnB1dDogdmVydGljZXMgb2YgYSBub24tZGVnZW5lcmF0ZSBzaW1w bGV4CgpmdW5jdGlvbiBbY2VudGVyLHJhZGl1c19zcXIsbGFtYmRhXSA9IG1pbmliYWxsX2Jm c3MgKHBvaW50c2V0KQoKIyMgQXV0aG9yOiBIYXJhbGQgQW5sYXVmIDxhbmxhdWYgYXQgaGVw IGRvdCB0dS1kYXJtc3RhZHQgZG90IGRlPgojIyBUaGlzIHNvZnR3YXJlIGlzIGRpc3RyaWJ1 dGVkIHVuZGVyIHRoZSB0ZXJtcyBvZiB0aGUgR1BMCgogIFtkLG1dID0gc2l6ZShwb2ludHNl dCk7CgogIGlmIChkID09IDAgfHwgbSA9PSAwKQogICAgZXJyb3IgKCdtaW5pYmFsbF9iZnNz OiBub25jb25mb3JtaW5nIGFyZ3VtZW50OiBwb2ludHNldCBtdXN0IGJlIGFycmF5IG9mIHNp emUgKGQsbSknKTsKICBlbmRpZgoKICBpZiAobSA+IGQrMSkKICAgIHdhcm5pbmcoJ21pbmli YWxsX2Jmc3M6IHRvbyBtYW55IHBvaW50cyBpbiBwb2ludHNldCwgdHJ1bmNhdGluZyBpbnB1 dC4uLicpOwogICAgbSA9IGQrMTsKICBlbmRpZgoKICAlJSBIYW5kbGUgdHJpdmlhbCBjYXNl cyBzZXBhcmF0ZWx5CiAgaWYgKG0gPT0gMSkKICAgIGNlbnRlciA9IHBvaW50c2V0KDosMSk7 CiAgICByYWRpdXNfc3FyID0gMDsKICAgIGxhbWJkYSA9IFsgMSBdOwogICAgcmV0dXJuCiAg ZWxzZWlmIChtID09IDIpCiAgICBjZW50ZXIgPSAocG9pbnRzZXQoOiwxKSArIHBvaW50c2V0 KDosMikpICogMC41OwogICAgcmFkaXVzX3NxciA9IG5vcm0ocG9pbnRzZXQoOiwxKSAtIHBv aW50c2V0KDosMikpXjIgKiAwLjI1OwogICAgbGFtYmRhID0gWyAwLjUgOyAwLjUgXTsKICAg IHJldHVybgogIGVuZGlmCgogICUlIFNoaWZ0IG9yaWdpbiB0byB0aGUgcG9zaXRpb24gb2Yg dGhlIGxhc3QgcG9pbnQKICBwMCA9IHBvaW50c2V0KDosbSk7CiAgZm9yIGkgPSAxOm0tMQog ICAgcSg6LGkpID0gcG9pbnRzZXQoOixpKSAtIHAwKDopOwogIGVuZGZvcgogIEEgPSBxJyAq IHE7CiAgZm9yIGkgPSAxOm0tMQogICAgeChpKSA9IDAuNSAqIEEoaSxpKTsKICBlbmRmb3IK ICBsYW1iZGEgPSBBIFwgeDsKICBsYW1iZGEobSkgPSAxIC0gc3VtKGxhbWJkYSgxOm0tMSkp OwoKICBjKDE6ZCkgPSAwOwogIGZvciBpID0gMTptLTEKICAgIGMgPSBjICsgbGFtYmRhKGkp ICogcSg6LGkpOwogIGVuZGZvcgogIGNlbnRlciA9IHAwICsgYzsKICByYWRpdXNfc3FyID0g YycgKiBjOwoKZW5kZnVuY3Rpb24KCiMjIExvY2FsIFZhcmlhYmxlczoKIyMgbW9kZTpvY3Rh dmUKIyMgcGFnZS1kZWxpbWl0ZXI6Il4jIAwiCiMjIEVuZDoK --NHgP2=-bQA=_U4Sw=?9SSKsrLnwKdFMoWqGDjv--