从表面坐标推导出流形的法向量
Deduce normal vector of manifold from surface coordinates
亲爱的 Whosebug 用户,
我是一名计算化学家,我有一道几何题。我有一堆定义分子表面的坐标,我想导出该表面的向外法向量。当我看它时,表面似乎近似于流形的属性,尽管坐标点不是使用该框架明确导出的。我还必须明确指出,在一般情况下,分子表面并不总是凸包,并且可以具有一定程度的凹度。他们没有的是不连续性,表面是光滑的结构。但由于我不知道如何处理这些数学规范,所以我试图为一般问题设计一个算法。
作为初步说明,对于表面的每个点,我可以确定最近原子的位置。因此,对于每个点,我也有这些可用的 xyz 坐标。该算法采用以下形式:
1. 计算每个可用点之间的距离矩阵(不可避免地会缩放到点数的平方,但对于我使用 numpy 的情况来说仍然合理)
2.提取每个点的两个最近邻
3.使用这个点的三元组生成以每个点为中心的两个向量
4.根据这两个向量的叉积得到法向量然后归一化
5.计算点与其底层原子之间的向量
6.如果这个向量与法向量的夹角小于90°,这个向量是向内指向的,因此它被它的相反方向所取代
这个完整过程的结果还算可以,但是当我使用 matplotlib 目视检查结果时,仍然有各种向量与表面有些平行。这是水分子的 matplotlib 结果:
这是用于比较的水分子表面(您可以在其中找到下面的原子)。忽略表面的颜色编码,它是由表面电荷进行颜色编码的,与目前的讨论无关。
这个surface是第三方软件获取的,我没法获取代码。我只能想象它,我无法访问其中用于最终渲染的平滑程序。
正如图片所示,表面非常光滑,因此我希望法向量能够解释这一点 "smoothness",但它们做得并不完美。我需要法向量来实际反映表面的平滑度,因为当前法向量描绘的表面粗糙度对我基于这些法向量的后续计算质量有显着影响。有谁知道在合理的计算时间内我可以做什么来解决这个问题?
这是一个可以重现我的第一个图形的工作代码:
import math
import matplotlib.pyplot as plt
import numpy as np
from mpl_toolkits.mplot3d import Axes3D
from sklearn.metrics.pairwise import euclidean_distances
coord = np.array([[ -1.2873729481345813 , 0.03256614731449952 , 1.5416924157851974 ],
[ 0.01652948394293976 , 1.163819319621563 , 1.6678642152662158 ],
[ 0.2794741299695759 , -0.5617278297487918 , 2.0029980474538105 ],
[ -0.6610946883884103 , -1.520269012229838 , 0.8599487353786675 ],
[ -1.0054760643749452 , 1.3940480132966795 , 0.33773540172063504 ],
[ 1.4883666878869983 , 0.43621600821755363 , 1.1450836953714032 ],
[ 1.093267571959524 , -1.3195317617899807 , 0.5499808804367919 ],
[ -0.044527698166979345 , -1.2654977063529413 , -0.7625313980846089 ],
[ 0.5831126343857715 , 1.5608391229347571 , -0.025329599172539626 ],
[ -1.0148260960520252 , 0.5209590347086124 , 1.6888058951547953 ],
[ -0.5081521680198725 , 0.9028613364971467 , 1.774471036317154 ],
[ -0.8515308971351676 , -0.19088994302497722 , 1.8836904989342724 ],
[ -0.28882376105739577 , -0.3548423909103548 , 2.05954224229103 ],
[ -1.2492850863979417 , -0.6364567978568106 , 1.3978142132864995 ],
[ -1.041946944653105 , -1.1796586713507826 , 1.095177216946184 ],
[ -1.5436272899575374 , -0.22919366151705667 , 1.1247655324014234 ],
[ -1.4689928122303186 , 0.5000569131417527 , 1.143407587573163 ],
[ -1.290329873679581 , 1.0646014214476844 , 0.8016157211074683 ],
[ 0.14732873180147785 , 0.6686847769257502 , 1.979342311827131 ],
[ 0.2131468214849169 , 0.056951984120299164 , 2.1073021163341092 ],
[ -0.0337308931325195 , -0.9942670995597854 , 1.8046165534409335 ],
[ -0.39947905888129415 , -1.37104542496434 , 1.3601926538530802 ],
[ -1.0634807007597644 , -1.3502061551644402 , 0.46773962277667314 ],
[ -0.7567553825092089 , 1.504835877060738 , 0.749651591341389 ],
[ -0.4472405474525535 , 1.4079282750475794 , 1.2824860152857014 ],
[ 1.5587883382335772 , -0.17395630632391745 , 1.1074452857884236 ],
[ 1.4235308457149791 , -0.8271279096055679 , 0.8993584919303867 ],
[ 0.7692195046037488 , -1.5167605398337778 , 0.1442407045573779 ],
[ 0.32445030395311525 , -1.4651179890494386 , -0.4390957189920336 ],
[ 0.00020957327211999695 , -0.9406835964416262 , -1.0384819480983047 ],
[ -0.025419507383719626 , 1.1326596798290633 , -0.892662493220727 ],
[ 0.273179306851356 , 1.4187473807855993 , -0.5317469722738123 ],
[ 1.0466981946916247 , 1.36444283710828 , 0.43534667343027367 ],
[ 1.3436060983763205 , 0.9793761802108056 , 0.8419277824771877 ],
[ 0.6151824124478109 , 0.9956899926113054 , 1.6618898508556756 ],
[ 1.1365292964079234 , 0.8011459463250883 , 1.4138732980146593 ],
[ 1.2547802693224217 , 0.08060290685487881 , 1.575156048146257 ],
[ 0.8061825016438882 , -0.24826596542943638 , 1.9004575307208522 ],
[ 0.7233645669036894 , -0.8826763855961071 , 1.6883835169082952 ],
[ 0.9497767017034661 , -1.2104420562255824 , 1.1703947277344628 ],
[ 0.5485586943697519 , -1.5967100787262565 , 0.7302017476622693 ],
[ -0.057205092501839166 , -1.6791781005999753 , 0.7696455114375087 ],
[ -0.48640523757119286 , -1.644761057333236 , 0.272618697792796 ],
[ -0.25748236328111623 , -1.501405645515218 , -0.39720000180351417 ],
[ -0.548591180200772 , 1.5966973503597166 , 0.07283122808381894 ],
[ 0.03137820171609954 , 1.6808890342631952 , 0.03811707406017944 ],
[ 0.5162256751875325 , 1.631621931751996 , 0.5739653241581716 ],
[ 0.29520943803145566 , 1.496670085663698 , 1.1960154986767626 ],
[ -0.4357510193390936 , 0.3045200764909755 , 2.0372933983710304 ],
[ -0.7252024085145093 , -0.8680072662231473 , 1.697293877794835 ],
[ -1.4026930025843594 , -0.8680775688445073 , 0.8886609086751069 ],
[ -1.0256665134561849 , 1.0589029487344046 , 1.2875897568848411 ],
[ 0.7343131121864293 , 0.4025541238612941 , 1.9038880538445522 ],
[ 0.2960177800611157 , -1.4085266358073394 , 1.3432377515245404 ],
[ -0.612740309720231 , -1.5270987389589377 , -0.09945502108597855 ],
[ -0.1766515084369774 , 1.6876887908829954 , 0.68242908228023 ],
[ 1.2576856987740817 , -0.6008401432185712 , 1.4092999027282396 ],
[ 0.20077831851211708 , -1.6825464196730353 , 0.10625820785053844 ],
[ -0.27284992140625597 , 1.4578679023663585 , -0.46947226287431315 ],
[ 0.9025274704849868 , 1.2861804462963813 , 1.101223178352364 ],
[ 1.7839350486036138 , -0.019389958495239716 , -1.0076276425199053 ],
[ 1.7377417775538346 , -0.8515596284340875 , -0.06551032812067904 ],
[ 1.9308978238593717 , 0.4526510335304934 , 0.15333098082293778 ],
[ 1.2640412923943216 , 1.1135989051613437 , -0.6487291165436305 ],
[ 0.7110361800083296 , 0.3011211057247756 , -1.4642623430903987 ],
[ 0.9695876990906258 , -0.9216814194628465 , -1.0943999970506841 ],
[ 2.04200929018949 , -0.17973742001943738 , -0.3635312878501747 ],
[ 1.8622593656466528 , 0.5635148250993117 , -0.6107377370213111 ],
[ 1.3955962627400398 , 0.5519618371910718 , -1.1944706185224625 ],
[ 1.2607629333920816 , -0.21551781723753682 , -1.3829616822422597 ],
[ 1.6583970787010358 , -0.6949731322511498 , -0.8399375243260477 ],
[ 1.9444143764920114 , -0.2454332569517364 , 0.2874392721860158 ],
[ 1.636157003293636 , 0.961650115510226 , -0.12326332301091819 ],
[ 0.8313542025004879 , 0.894165587734687 , -1.1420580452198033 ],
[ 0.6430205539552906 , -0.46024492898875324 , -1.4104511496936394 ],
[ 1.2947193164727608 , -1.1400972481068032 , -0.5315305572334722 ],
[ -1.7839430005914738 , 0.019376416779039715 , -1.0076136023161453 ],
[ -1.7377420093346745 , 0.8515508599214876 , -0.06550088225767904 ],
[ -1.9308962617200118 , -0.45265869341099335 , 0.15334861627561777 ],
[ -1.2640463026705615 , -1.1136106280858835 , -0.6487135972817705 ],
[ -0.7110478738279696 , -0.3011369604867556 , -1.4642554706297384 ],
[ -0.9695963617672257 , 0.9216674385272464 , -1.0943972008635638 ],
[ -2.04201196413603 , 0.17972714175629736 , -0.36351594480525473 ],
[ -1.8622640647650528 , -0.5635263554023318 , -0.6107201020978111 ],
[ -1.3956057456456397 , -0.5519763250811119 , -1.1944568661926225 ],
[ -1.2607739609741015 , 0.21550237470677686 , -1.3829529227257198 ],
[ -1.6584036564084357 , 0.6949604403980297 , -0.8399279355844478 ],
[ -1.9444117157749716 , 0.24542627653835639 , 0.2874534822565558 ],
[ -1.636157708161396 , -0.961659176659366 , -0.12324552404161819 ],
[ -0.8313632557119278 , -0.8941798105055468 , -1.1420471832711232 ],
[ -0.6430318069679907 , 0.46022934675447325 , -1.4104486921817194 ],
[ -1.294723366816481 , 1.1400861183930433 , -0.5315262031404322 ],
[ -4.5154929399999335e-06 , -0.6700687207417302 , -1.1844736267236027 ],
[ 0.16109703335223763 , -0.6271001397877708 , -1.2186543868869022 ],
[ -0.16110646810245766 , -0.6271000117262108 , -1.2186531586601221 ],
[ -4.0048342399999414e-06 , 0.6700542836529702 , -1.1844770214133027 ],
[ 0.16109751120177762 , 0.6270854068873908 , -1.218657563554442 ],
[ -0.16110599025291766 , 0.6270855243653509 , -1.2186563353276623 ],
[ 0.36994498278851456 , 0.7707305304675487 , -1.169504657558063 ],
[ 0.7025740208231097 , 1.1378228594549835 , -0.7470950439135691 ],
[ 1.1902461772283626 , 1.1877239069058625 , -0.12779283816255813 ],
[ 1.5961784560313765 , 0.748033999209189 , 0.38770812025235435 ],
[ 1.7530425415210142 , -3.3999814999999503e-06 , 0.5869144868197315 ],
[ 1.5961778861045166 , -0.748041688194589 , 0.3877119097103343 ],
[ 1.1902452718013825 , -1.1877338978242624 , -0.12778682138595815 ],
[ 0.7025731540262697 , -1.1378356163972434 , -0.7470892795558292 ],
[ 0.3699443953987146 , -0.7707451739365088 , -1.1695007532680228 ],
[ 0.27881228877701597 , 1.0004381797559854 , -0.9277686907765464 ],
[ 0.7663313983885688 , 1.3243776625508408 , -0.3086616683059155 ],
[ 1.2752305145620413 , 1.102396662519724 , 0.337597800103595 ],
[ 1.5957376215740167 , 0.42599740483075377 , 0.7446167357487491 ],
[ 1.5957372971866766 , -0.4260032850789138 , 0.7446188937447891 ],
[ 1.2752296742242015 , -1.102404360501184 , 0.3376033845401351 ],
[ 0.7663303887131288 , -1.3243882472092006 , -0.3086549593618755 ],
[ 0.27881152675781595 , -1.0004515288506655 , -0.9277636228196864 ],
[ -0.36995352216617455 , -0.7707449252219087 , -1.169497906808803 ],
[ -0.7025791448730497 , -1.1378352248040433 , -0.747083649080629 ],
[ -1.1902463164027026 , -1.1877331892522427 , -0.12777707971133812 ],
[ -1.5961744570181167 , -0.748040609725749 , 0.3877250868215143 ],
[ -1.7530369449133343 , -2.0638019999999696e-06 , 0.5869289937602514 ],
[ -1.5961738870912567 , 0.748035353909989 , 0.38772129736353433 ],
[ -1.1902454109757226 , 1.1877250123628826 , -0.12778309701711812 ],
[ -0.7025782780762097 , 1.1378235389221032 , -0.747089413438369 ],
[ -0.3699529347763746 , 0.7707308453296488 , -1.1695018110988429 ],
[ -0.27881865163733593 , -1.0004514007891054 , -0.9277613558125664 ],
[ -0.7663327657896889 , -1.3243878630245205 , -0.3086485986182755 ],
[ -1.2752266880614613 , -1.102403584723304 , 0.33761404540041506 ],
[ -1.5957305300328366 , -0.42600214945863374 , 0.7446322698276491 ],
[ -1.5957302056454967 , 0.42599870132175377 , 0.7446301118316091 ],
[ -1.2752258482528014 , 1.102397829890804 , 0.33760846043469506 ],
[ -0.7663317561142488 , 1.3243784467956008 , -0.3086553075623155 ],
[ -0.2788178890889559 , 1.0004384766259653 , -0.9277664242986065 ],
[ -0.23018790924333662 , 0.4635384934491132 , -1.3322307085678806 ],
[ -0.23018827914015663 , -0.20734291124417695 , -1.3622983801385802 ],
[ 0.23017815010577664 , 0.20732757296187695 , -1.3623011720922602 ],
[ 0.23017800828553664 , -0.46355367932757324 , -1.3322301015984204 ]])
atoms = np.asarray([[ -1.21073950301e-06 , 0.0 , 0.40150605663 ],
[ -1.21073950301e-06 , 0.0 , 0.40150605663 ],
[ -1.21073950301e-06 , 0.0 , 0.40150605663 ],
[ -1.21073950301e-06 , 0.0 , 0.40150605663 ],
[ -1.21073950301e-06 , 0.0 , 0.40150605663 ],
[ -1.21073950301e-06 , 0.0 , 0.40150605663 ],
[ -1.21073950301e-06 , 0.0 , 0.40150605663 ],
[ -1.21073950301e-06 , 0.0 , 0.40150605663 ],
[ -1.21073950301e-06 , 0.0 , 0.40150605663 ],
[ -1.21073950301e-06 , 0.0 , 0.40150605663 ],
[ -1.21073950301e-06 , 0.0 , 0.40150605663 ],
[ -1.21073950301e-06 , 0.0 , 0.40150605663 ],
[ -1.21073950301e-06 , 0.0 , 0.40150605663 ],
[ -1.21073950301e-06 , 0.0 , 0.40150605663 ],
[ -1.21073950301e-06 , 0.0 , 0.40150605663 ],
[ -1.21073950301e-06 , 0.0 , 0.40150605663 ],
[ -1.21073950301e-06 , 0.0 , 0.40150605663 ],
[ -1.21073950301e-06 , 0.0 , 0.40150605663 ],
[ -1.21073950301e-06 , 0.0 , 0.40150605663 ],
[ -1.21073950301e-06 , 0.0 , 0.40150605663 ],
[ -1.21073950301e-06 , 0.0 , 0.40150605663 ],
[ -1.21073950301e-06 , 0.0 , 0.40150605663 ],
[ -1.21073950301e-06 , 0.0 , 0.40150605663 ],
[ -1.21073950301e-06 , 0.0 , 0.40150605663 ],
[ -1.21073950301e-06 , 0.0 , 0.40150605663 ],
[ -1.21073950301e-06 , 0.0 , 0.40150605663 ],
[ -1.21073950301e-06 , 0.0 , 0.40150605663 ],
[ -1.21073950301e-06 , 0.0 , 0.40150605663 ],
[ -1.21073950301e-06 , 0.0 , 0.40150605663 ],
[ -1.21073950301e-06 , 0.0 , 0.40150605663 ],
[ -1.21073950301e-06 , 0.0 , 0.40150605663 ],
[ -1.21073950301e-06 , 0.0 , 0.40150605663 ],
[ -1.21073950301e-06 , 0.0 , 0.40150605663 ],
[ -1.21073950301e-06 , 0.0 , 0.40150605663 ],
[ -1.21073950301e-06 , 0.0 , 0.40150605663 ],
[ -1.21073950301e-06 , 0.0 , 0.40150605663 ],
[ -1.21073950301e-06 , 0.0 , 0.40150605663 ],
[ -1.21073950301e-06 , 0.0 , 0.40150605663 ],
[ -1.21073950301e-06 , 0.0 , 0.40150605663 ],
[ -1.21073950301e-06 , 0.0 , 0.40150605663 ],
[ -1.21073950301e-06 , 0.0 , 0.40150605663 ],
[ -1.21073950301e-06 , 0.0 , 0.40150605663 ],
[ -1.21073950301e-06 , 0.0 , 0.40150605663 ],
[ -1.21073950301e-06 , 0.0 , 0.40150605663 ],
[ -1.21073950301e-06 , 0.0 , 0.40150605663 ],
[ -1.21073950301e-06 , 0.0 , 0.40150605663 ],
[ -1.21073950301e-06 , 0.0 , 0.40150605663 ],
[ -1.21073950301e-06 , 0.0 , 0.40150605663 ],
[ -1.21073950301e-06 , 0.0 , 0.40150605663 ],
[ -1.21073950301e-06 , 0.0 , 0.40150605663 ],
[ -1.21073950301e-06 , 0.0 , 0.40150605663 ],
[ -1.21073950301e-06 , 0.0 , 0.40150605663 ],
[ -1.21073950301e-06 , 0.0 , 0.40150605663 ],
[ -1.21073950301e-06 , 0.0 , 0.40150605663 ],
[ -1.21073950301e-06 , 0.0 , 0.40150605663 ],
[ -1.21073950301e-06 , 0.0 , 0.40150605663 ],
[ -1.21073950301e-06 , 0.0 , 0.40150605663 ],
[ -1.21073950301e-06 , 0.0 , 0.40150605663 ],
[ -1.21073950301e-06 , 0.0 , 0.40150605663 ],
[ -1.21073950301e-06 , 0.0 , 0.40150605663 ],
[ 0.764816906209 , 0.0 , -0.200753505843 ],
[ 0.764816906209 , 0.0 , -0.200753505843 ],
[ 0.764816906209 , 0.0 , -0.200753505843 ],
[ 0.764816906209 , 0.0 , -0.200753505843 ],
[ 0.764816906209 , 0.0 , -0.200753505843 ],
[ 0.764816906209 , 0.0 , -0.200753505843 ],
[ 0.764816906209 , 0.0 , -0.200753505843 ],
[ 0.764816906209 , 0.0 , -0.200753505843 ],
[ 0.764816906209 , 0.0 , -0.200753505843 ],
[ 0.764816906209 , 0.0 , -0.200753505843 ],
[ 0.764816906209 , 0.0 , -0.200753505843 ],
[ 0.764816906209 , 0.0 , -0.200753505843 ],
[ 0.764816906209 , 0.0 , -0.200753505843 ],
[ 0.764816906209 , 0.0 , -0.200753505843 ],
[ 0.764816906209 , 0.0 , -0.200753505843 ],
[ 0.764816906209 , 0.0 , -0.200753505843 ],
[ -0.764815694412 , 0.0 , -0.200752550787 ],
[ -0.764815694412 , 0.0 , -0.200752550787 ],
[ -0.764815694412 , 0.0 , -0.200752550787 ],
[ -0.764815694412 , 0.0 , -0.200752550787 ],
[ -0.764815694412 , 0.0 , -0.200752550787 ],
[ -0.764815694412 , 0.0 , -0.200752550787 ],
[ -0.764815694412 , 0.0 , -0.200752550787 ],
[ -0.764815694412 , 0.0 , -0.200752550787 ],
[ -0.764815694412 , 0.0 , -0.200752550787 ],
[ -0.764815694412 , 0.0 , -0.200752550787 ],
[ -0.764815694412 , 0.0 , -0.200752550787 ],
[ -0.764815694412 , 0.0 , -0.200752550787 ],
[ -0.764815694412 , 0.0 , -0.200752550787 ],
[ -0.764815694412 , 0.0 , -0.200752550787 ],
[ -0.764815694412 , 0.0 , -0.200752550787 ],
[ -0.764815694412 , 0.0 , -0.200752550787 ],
[ -1.21073950301e-06 , 0.0 , 0.40150605663 ],
[ 0.764816906209 , 0.0 , -0.200753505843 ],
[ -0.764815694412 , 0.0 , -0.200752550787 ],
[ -1.21073950301e-06 , 0.0 , 0.40150605663 ],
[ 0.764816906209 , 0.0 , -0.200753505843 ],
[ -0.764815694412 , 0.0 , -0.200752550787 ],
[ 0.764816906209 , 0.0 , -0.200753505843 ],
[ 0.764816906209 , 0.0 , -0.200753505843 ],
[ 0.764816906209 , 0.0 , -0.200753505843 ],
[ 0.764816906209 , 0.0 , -0.200753505843 ],
[ 0.764816906209 , 0.0 , -0.200753505843 ],
[ 0.764816906209 , 0.0 , -0.200753505843 ],
[ 0.764816906209 , 0.0 , -0.200753505843 ],
[ 0.764816906209 , 0.0 , -0.200753505843 ],
[ 0.764816906209 , 0.0 , -0.200753505843 ],
[ -1.21073950301e-06 , 0.0 , 0.40150605663 ],
[ -1.21073950301e-06 , 0.0 , 0.40150605663 ],
[ -1.21073950301e-06 , 0.0 , 0.40150605663 ],
[ -1.21073950301e-06 , 0.0 , 0.40150605663 ],
[ -1.21073950301e-06 , 0.0 , 0.40150605663 ],
[ -1.21073950301e-06 , 0.0 , 0.40150605663 ],
[ -1.21073950301e-06 , 0.0 , 0.40150605663 ],
[ -1.21073950301e-06 , 0.0 , 0.40150605663 ],
[ -0.764815694412 , 0.0 , -0.200752550787 ],
[ -0.764815694412 , 0.0 , -0.200752550787 ],
[ -0.764815694412 , 0.0 , -0.200752550787 ],
[ -0.764815694412 , 0.0 , -0.200752550787 ],
[ -0.764815694412 , 0.0 , -0.200752550787 ],
[ -0.764815694412 , 0.0 , -0.200752550787 ],
[ -0.764815694412 , 0.0 , -0.200752550787 ],
[ -0.764815694412 , 0.0 , -0.200752550787 ],
[ -0.764815694412 , 0.0 , -0.200752550787 ],
[ -1.21073950301e-06 , 0.0 , 0.40150605663 ],
[ -1.21073950301e-06 , 0.0 , 0.40150605663 ],
[ -1.21073950301e-06 , 0.0 , 0.40150605663 ],
[ -1.21073950301e-06 , 0.0 , 0.40150605663 ],
[ -1.21073950301e-06 , 0.0 , 0.40150605663 ],
[ -1.21073950301e-06 , 0.0 , 0.40150605663 ],
[ -1.21073950301e-06 , 0.0 , 0.40150605663 ],
[ -1.21073950301e-06 , 0.0 , 0.40150605663 ],
[ -0.764815694412 , 0.0 , -0.200752550787 ],
[ -0.764815694412 , 0.0 , -0.200752550787 ],
[ 0.764816906209 , 0.0 , -0.200753505843 ],
[ 0.764816906209 , 0.0 , -0.200753505843 ]])
#transpose coordinate array
xyz = np.transpose(coord)
#establish distance matrix
n = len(coord)
dist_vectors = list()
dist_vectors.append(xyz[0])
dist_vectors.append(xyz[1])
dist_vectors.append(xyz[2])
dist_vectors = map(list, zip(*dist_vectors))
d = euclidean_distances(dist_vectors, dist_vectors)
#find two nearest neighbors for each segment
x1 = np.zeros((n))
x2 = np.zeros((n))
y1 = np.zeros((n))
y2 = np.zeros((n))
z1 = np.zeros((n))
z2 = np.zeros((n))
for i in range(n):
x_copy = xyz[0]
y_copy = xyz[1]
z_copy = xyz[2]
d1 = np.delete(d[i], i) #removes distance between segment and itself
x_copy = np.delete(x_copy, i)
y_copy = np.delete(y_copy, i)
z_copy = np.delete(z_copy, i)
j1 = np.argmin(d1) #get indice of minimum distance
x1[i] = x_copy[j1]
y1[i] = y_copy[j1]
z1[i] = z_copy[j1]
d2 = np.delete(d1, j1) #removes minimum distance
x_copy = np.delete(x_copy, j1)
y_copy = np.delete(y_copy, j1)
z_copy = np.delete(z_copy, j1)
j2 = np.argmin(d2) #get indice of second minimum distance
x2[i] = x_copy[j2]
y2[i] = y_copy[j2]
z2[i] = z_copy[j2]
#compute normal vector for each segment based on cross product
normal = list()
forGraphs = list()
for i in range(n):
#make vectors for cross product
v1 = np.zeros((3))
v1[0] = x1[i] - coord[i][0]
v1[1] = y1[i] - coord[i][1]
v1[2] = z1[i] - coord[i][2]
v2 = np.zeros((3))
v2[0] = x2[i] - coord[i][0]
v2[1] = y2[i] - coord[i][1]
v2[2] = z2[i] - coord[i][2]
#make cross product and normalize (normal vector should have a unit norm)
nv = np.cross(v1, v2)
nv = nv / np.linalg.norm(nv)
normal.append(nv)
#check of outwards pointing
atv = np.zeros((3))
atv[0] = atoms[i][0] - coord[i][0]
atv[1] = atoms[i][1] - coord[i][1]
atv[2] = atoms[i][2] - coord[i][2]
th_check = math.acos(np.dot(nv, atv) / (np.linalg.norm(nv) * np.linalg.norm(atv)))
if th_check < (math.pi / 2): #if inwards pointing (i. e. pointing towards underlying atom), normal vector is replaced by its opposite
nv[0] = -nv[0]
nv[1] = -nv[1]
nv[2] = -nv[2]
forGraphs.append(np.array([coord[i][0],coord[i][1],coord[i][2],nv[0],nv[1], nv[2]]))
#plot normal vectors (for checkup)
forGraphs = np.asarray(forGraphs)
X, Y, Z, U, V, W = zip(*forGraphs)
fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
ax.quiver(X, Y, Z, U, V, W)
ax.set_xlim([min(xyz[0])- 1, max(xyz[0]) + 1])
ax.set_ylim([min(xyz[1])- 1, max(xyz[1]) + 1])
ax.set_zlim([min(xyz[2])- 1, max(xyz[2]) + 1])
plt.show()
此代码的前 280 行主要用于重现结果所需的坐标表。这段代码最重要的部分是第 282 行到第 355 行,其中实现了我刚刚概述的算法。
在此先感谢您的帮助!
估计法线的"standard"过程是为每个点找到k个最近的邻居,其中k是一个小数字(十?)。然后计算通过这些点的最佳平面拟合,并使用平面的法线。
不幸的是,法线的符号是不确定的,你需要实施一个法线一致性强制过程。也许在你的情况下这更容易,因为所有法线似乎都指向远离某个中心点。
我按照@Yves Daoust 的建议使用了最近邻技术,并使用奇异值分解 (SVD) 通过 10 个最近邻来拟合一个平面。在视觉上简单的情况下,这仍然会给出错误的答案。我不明白为什么。这是我得到的示例:
亲爱的 Whosebug 用户,
我是一名计算化学家,我有一道几何题。我有一堆定义分子表面的坐标,我想导出该表面的向外法向量。当我看它时,表面似乎近似于流形的属性,尽管坐标点不是使用该框架明确导出的。我还必须明确指出,在一般情况下,分子表面并不总是凸包,并且可以具有一定程度的凹度。他们没有的是不连续性,表面是光滑的结构。但由于我不知道如何处理这些数学规范,所以我试图为一般问题设计一个算法。
作为初步说明,对于表面的每个点,我可以确定最近原子的位置。因此,对于每个点,我也有这些可用的 xyz 坐标。该算法采用以下形式:
1. 计算每个可用点之间的距离矩阵(不可避免地会缩放到点数的平方,但对于我使用 numpy 的情况来说仍然合理)
2.提取每个点的两个最近邻
3.使用这个点的三元组生成以每个点为中心的两个向量
4.根据这两个向量的叉积得到法向量然后归一化
5.计算点与其底层原子之间的向量
6.如果这个向量与法向量的夹角小于90°,这个向量是向内指向的,因此它被它的相反方向所取代
这个完整过程的结果还算可以,但是当我使用 matplotlib 目视检查结果时,仍然有各种向量与表面有些平行。这是水分子的 matplotlib 结果:
这个surface是第三方软件获取的,我没法获取代码。我只能想象它,我无法访问其中用于最终渲染的平滑程序。
正如图片所示,表面非常光滑,因此我希望法向量能够解释这一点 "smoothness",但它们做得并不完美。我需要法向量来实际反映表面的平滑度,因为当前法向量描绘的表面粗糙度对我基于这些法向量的后续计算质量有显着影响。有谁知道在合理的计算时间内我可以做什么来解决这个问题?
这是一个可以重现我的第一个图形的工作代码:
import math
import matplotlib.pyplot as plt
import numpy as np
from mpl_toolkits.mplot3d import Axes3D
from sklearn.metrics.pairwise import euclidean_distances
coord = np.array([[ -1.2873729481345813 , 0.03256614731449952 , 1.5416924157851974 ],
[ 0.01652948394293976 , 1.163819319621563 , 1.6678642152662158 ],
[ 0.2794741299695759 , -0.5617278297487918 , 2.0029980474538105 ],
[ -0.6610946883884103 , -1.520269012229838 , 0.8599487353786675 ],
[ -1.0054760643749452 , 1.3940480132966795 , 0.33773540172063504 ],
[ 1.4883666878869983 , 0.43621600821755363 , 1.1450836953714032 ],
[ 1.093267571959524 , -1.3195317617899807 , 0.5499808804367919 ],
[ -0.044527698166979345 , -1.2654977063529413 , -0.7625313980846089 ],
[ 0.5831126343857715 , 1.5608391229347571 , -0.025329599172539626 ],
[ -1.0148260960520252 , 0.5209590347086124 , 1.6888058951547953 ],
[ -0.5081521680198725 , 0.9028613364971467 , 1.774471036317154 ],
[ -0.8515308971351676 , -0.19088994302497722 , 1.8836904989342724 ],
[ -0.28882376105739577 , -0.3548423909103548 , 2.05954224229103 ],
[ -1.2492850863979417 , -0.6364567978568106 , 1.3978142132864995 ],
[ -1.041946944653105 , -1.1796586713507826 , 1.095177216946184 ],
[ -1.5436272899575374 , -0.22919366151705667 , 1.1247655324014234 ],
[ -1.4689928122303186 , 0.5000569131417527 , 1.143407587573163 ],
[ -1.290329873679581 , 1.0646014214476844 , 0.8016157211074683 ],
[ 0.14732873180147785 , 0.6686847769257502 , 1.979342311827131 ],
[ 0.2131468214849169 , 0.056951984120299164 , 2.1073021163341092 ],
[ -0.0337308931325195 , -0.9942670995597854 , 1.8046165534409335 ],
[ -0.39947905888129415 , -1.37104542496434 , 1.3601926538530802 ],
[ -1.0634807007597644 , -1.3502061551644402 , 0.46773962277667314 ],
[ -0.7567553825092089 , 1.504835877060738 , 0.749651591341389 ],
[ -0.4472405474525535 , 1.4079282750475794 , 1.2824860152857014 ],
[ 1.5587883382335772 , -0.17395630632391745 , 1.1074452857884236 ],
[ 1.4235308457149791 , -0.8271279096055679 , 0.8993584919303867 ],
[ 0.7692195046037488 , -1.5167605398337778 , 0.1442407045573779 ],
[ 0.32445030395311525 , -1.4651179890494386 , -0.4390957189920336 ],
[ 0.00020957327211999695 , -0.9406835964416262 , -1.0384819480983047 ],
[ -0.025419507383719626 , 1.1326596798290633 , -0.892662493220727 ],
[ 0.273179306851356 , 1.4187473807855993 , -0.5317469722738123 ],
[ 1.0466981946916247 , 1.36444283710828 , 0.43534667343027367 ],
[ 1.3436060983763205 , 0.9793761802108056 , 0.8419277824771877 ],
[ 0.6151824124478109 , 0.9956899926113054 , 1.6618898508556756 ],
[ 1.1365292964079234 , 0.8011459463250883 , 1.4138732980146593 ],
[ 1.2547802693224217 , 0.08060290685487881 , 1.575156048146257 ],
[ 0.8061825016438882 , -0.24826596542943638 , 1.9004575307208522 ],
[ 0.7233645669036894 , -0.8826763855961071 , 1.6883835169082952 ],
[ 0.9497767017034661 , -1.2104420562255824 , 1.1703947277344628 ],
[ 0.5485586943697519 , -1.5967100787262565 , 0.7302017476622693 ],
[ -0.057205092501839166 , -1.6791781005999753 , 0.7696455114375087 ],
[ -0.48640523757119286 , -1.644761057333236 , 0.272618697792796 ],
[ -0.25748236328111623 , -1.501405645515218 , -0.39720000180351417 ],
[ -0.548591180200772 , 1.5966973503597166 , 0.07283122808381894 ],
[ 0.03137820171609954 , 1.6808890342631952 , 0.03811707406017944 ],
[ 0.5162256751875325 , 1.631621931751996 , 0.5739653241581716 ],
[ 0.29520943803145566 , 1.496670085663698 , 1.1960154986767626 ],
[ -0.4357510193390936 , 0.3045200764909755 , 2.0372933983710304 ],
[ -0.7252024085145093 , -0.8680072662231473 , 1.697293877794835 ],
[ -1.4026930025843594 , -0.8680775688445073 , 0.8886609086751069 ],
[ -1.0256665134561849 , 1.0589029487344046 , 1.2875897568848411 ],
[ 0.7343131121864293 , 0.4025541238612941 , 1.9038880538445522 ],
[ 0.2960177800611157 , -1.4085266358073394 , 1.3432377515245404 ],
[ -0.612740309720231 , -1.5270987389589377 , -0.09945502108597855 ],
[ -0.1766515084369774 , 1.6876887908829954 , 0.68242908228023 ],
[ 1.2576856987740817 , -0.6008401432185712 , 1.4092999027282396 ],
[ 0.20077831851211708 , -1.6825464196730353 , 0.10625820785053844 ],
[ -0.27284992140625597 , 1.4578679023663585 , -0.46947226287431315 ],
[ 0.9025274704849868 , 1.2861804462963813 , 1.101223178352364 ],
[ 1.7839350486036138 , -0.019389958495239716 , -1.0076276425199053 ],
[ 1.7377417775538346 , -0.8515596284340875 , -0.06551032812067904 ],
[ 1.9308978238593717 , 0.4526510335304934 , 0.15333098082293778 ],
[ 1.2640412923943216 , 1.1135989051613437 , -0.6487291165436305 ],
[ 0.7110361800083296 , 0.3011211057247756 , -1.4642623430903987 ],
[ 0.9695876990906258 , -0.9216814194628465 , -1.0943999970506841 ],
[ 2.04200929018949 , -0.17973742001943738 , -0.3635312878501747 ],
[ 1.8622593656466528 , 0.5635148250993117 , -0.6107377370213111 ],
[ 1.3955962627400398 , 0.5519618371910718 , -1.1944706185224625 ],
[ 1.2607629333920816 , -0.21551781723753682 , -1.3829616822422597 ],
[ 1.6583970787010358 , -0.6949731322511498 , -0.8399375243260477 ],
[ 1.9444143764920114 , -0.2454332569517364 , 0.2874392721860158 ],
[ 1.636157003293636 , 0.961650115510226 , -0.12326332301091819 ],
[ 0.8313542025004879 , 0.894165587734687 , -1.1420580452198033 ],
[ 0.6430205539552906 , -0.46024492898875324 , -1.4104511496936394 ],
[ 1.2947193164727608 , -1.1400972481068032 , -0.5315305572334722 ],
[ -1.7839430005914738 , 0.019376416779039715 , -1.0076136023161453 ],
[ -1.7377420093346745 , 0.8515508599214876 , -0.06550088225767904 ],
[ -1.9308962617200118 , -0.45265869341099335 , 0.15334861627561777 ],
[ -1.2640463026705615 , -1.1136106280858835 , -0.6487135972817705 ],
[ -0.7110478738279696 , -0.3011369604867556 , -1.4642554706297384 ],
[ -0.9695963617672257 , 0.9216674385272464 , -1.0943972008635638 ],
[ -2.04201196413603 , 0.17972714175629736 , -0.36351594480525473 ],
[ -1.8622640647650528 , -0.5635263554023318 , -0.6107201020978111 ],
[ -1.3956057456456397 , -0.5519763250811119 , -1.1944568661926225 ],
[ -1.2607739609741015 , 0.21550237470677686 , -1.3829529227257198 ],
[ -1.6584036564084357 , 0.6949604403980297 , -0.8399279355844478 ],
[ -1.9444117157749716 , 0.24542627653835639 , 0.2874534822565558 ],
[ -1.636157708161396 , -0.961659176659366 , -0.12324552404161819 ],
[ -0.8313632557119278 , -0.8941798105055468 , -1.1420471832711232 ],
[ -0.6430318069679907 , 0.46022934675447325 , -1.4104486921817194 ],
[ -1.294723366816481 , 1.1400861183930433 , -0.5315262031404322 ],
[ -4.5154929399999335e-06 , -0.6700687207417302 , -1.1844736267236027 ],
[ 0.16109703335223763 , -0.6271001397877708 , -1.2186543868869022 ],
[ -0.16110646810245766 , -0.6271000117262108 , -1.2186531586601221 ],
[ -4.0048342399999414e-06 , 0.6700542836529702 , -1.1844770214133027 ],
[ 0.16109751120177762 , 0.6270854068873908 , -1.218657563554442 ],
[ -0.16110599025291766 , 0.6270855243653509 , -1.2186563353276623 ],
[ 0.36994498278851456 , 0.7707305304675487 , -1.169504657558063 ],
[ 0.7025740208231097 , 1.1378228594549835 , -0.7470950439135691 ],
[ 1.1902461772283626 , 1.1877239069058625 , -0.12779283816255813 ],
[ 1.5961784560313765 , 0.748033999209189 , 0.38770812025235435 ],
[ 1.7530425415210142 , -3.3999814999999503e-06 , 0.5869144868197315 ],
[ 1.5961778861045166 , -0.748041688194589 , 0.3877119097103343 ],
[ 1.1902452718013825 , -1.1877338978242624 , -0.12778682138595815 ],
[ 0.7025731540262697 , -1.1378356163972434 , -0.7470892795558292 ],
[ 0.3699443953987146 , -0.7707451739365088 , -1.1695007532680228 ],
[ 0.27881228877701597 , 1.0004381797559854 , -0.9277686907765464 ],
[ 0.7663313983885688 , 1.3243776625508408 , -0.3086616683059155 ],
[ 1.2752305145620413 , 1.102396662519724 , 0.337597800103595 ],
[ 1.5957376215740167 , 0.42599740483075377 , 0.7446167357487491 ],
[ 1.5957372971866766 , -0.4260032850789138 , 0.7446188937447891 ],
[ 1.2752296742242015 , -1.102404360501184 , 0.3376033845401351 ],
[ 0.7663303887131288 , -1.3243882472092006 , -0.3086549593618755 ],
[ 0.27881152675781595 , -1.0004515288506655 , -0.9277636228196864 ],
[ -0.36995352216617455 , -0.7707449252219087 , -1.169497906808803 ],
[ -0.7025791448730497 , -1.1378352248040433 , -0.747083649080629 ],
[ -1.1902463164027026 , -1.1877331892522427 , -0.12777707971133812 ],
[ -1.5961744570181167 , -0.748040609725749 , 0.3877250868215143 ],
[ -1.7530369449133343 , -2.0638019999999696e-06 , 0.5869289937602514 ],
[ -1.5961738870912567 , 0.748035353909989 , 0.38772129736353433 ],
[ -1.1902454109757226 , 1.1877250123628826 , -0.12778309701711812 ],
[ -0.7025782780762097 , 1.1378235389221032 , -0.747089413438369 ],
[ -0.3699529347763746 , 0.7707308453296488 , -1.1695018110988429 ],
[ -0.27881865163733593 , -1.0004514007891054 , -0.9277613558125664 ],
[ -0.7663327657896889 , -1.3243878630245205 , -0.3086485986182755 ],
[ -1.2752266880614613 , -1.102403584723304 , 0.33761404540041506 ],
[ -1.5957305300328366 , -0.42600214945863374 , 0.7446322698276491 ],
[ -1.5957302056454967 , 0.42599870132175377 , 0.7446301118316091 ],
[ -1.2752258482528014 , 1.102397829890804 , 0.33760846043469506 ],
[ -0.7663317561142488 , 1.3243784467956008 , -0.3086553075623155 ],
[ -0.2788178890889559 , 1.0004384766259653 , -0.9277664242986065 ],
[ -0.23018790924333662 , 0.4635384934491132 , -1.3322307085678806 ],
[ -0.23018827914015663 , -0.20734291124417695 , -1.3622983801385802 ],
[ 0.23017815010577664 , 0.20732757296187695 , -1.3623011720922602 ],
[ 0.23017800828553664 , -0.46355367932757324 , -1.3322301015984204 ]])
atoms = np.asarray([[ -1.21073950301e-06 , 0.0 , 0.40150605663 ],
[ -1.21073950301e-06 , 0.0 , 0.40150605663 ],
[ -1.21073950301e-06 , 0.0 , 0.40150605663 ],
[ -1.21073950301e-06 , 0.0 , 0.40150605663 ],
[ -1.21073950301e-06 , 0.0 , 0.40150605663 ],
[ -1.21073950301e-06 , 0.0 , 0.40150605663 ],
[ -1.21073950301e-06 , 0.0 , 0.40150605663 ],
[ -1.21073950301e-06 , 0.0 , 0.40150605663 ],
[ -1.21073950301e-06 , 0.0 , 0.40150605663 ],
[ -1.21073950301e-06 , 0.0 , 0.40150605663 ],
[ -1.21073950301e-06 , 0.0 , 0.40150605663 ],
[ -1.21073950301e-06 , 0.0 , 0.40150605663 ],
[ -1.21073950301e-06 , 0.0 , 0.40150605663 ],
[ -1.21073950301e-06 , 0.0 , 0.40150605663 ],
[ -1.21073950301e-06 , 0.0 , 0.40150605663 ],
[ -1.21073950301e-06 , 0.0 , 0.40150605663 ],
[ -1.21073950301e-06 , 0.0 , 0.40150605663 ],
[ -1.21073950301e-06 , 0.0 , 0.40150605663 ],
[ -1.21073950301e-06 , 0.0 , 0.40150605663 ],
[ -1.21073950301e-06 , 0.0 , 0.40150605663 ],
[ -1.21073950301e-06 , 0.0 , 0.40150605663 ],
[ -1.21073950301e-06 , 0.0 , 0.40150605663 ],
[ -1.21073950301e-06 , 0.0 , 0.40150605663 ],
[ -1.21073950301e-06 , 0.0 , 0.40150605663 ],
[ -1.21073950301e-06 , 0.0 , 0.40150605663 ],
[ -1.21073950301e-06 , 0.0 , 0.40150605663 ],
[ -1.21073950301e-06 , 0.0 , 0.40150605663 ],
[ -1.21073950301e-06 , 0.0 , 0.40150605663 ],
[ -1.21073950301e-06 , 0.0 , 0.40150605663 ],
[ -1.21073950301e-06 , 0.0 , 0.40150605663 ],
[ -1.21073950301e-06 , 0.0 , 0.40150605663 ],
[ -1.21073950301e-06 , 0.0 , 0.40150605663 ],
[ -1.21073950301e-06 , 0.0 , 0.40150605663 ],
[ -1.21073950301e-06 , 0.0 , 0.40150605663 ],
[ -1.21073950301e-06 , 0.0 , 0.40150605663 ],
[ -1.21073950301e-06 , 0.0 , 0.40150605663 ],
[ -1.21073950301e-06 , 0.0 , 0.40150605663 ],
[ -1.21073950301e-06 , 0.0 , 0.40150605663 ],
[ -1.21073950301e-06 , 0.0 , 0.40150605663 ],
[ -1.21073950301e-06 , 0.0 , 0.40150605663 ],
[ -1.21073950301e-06 , 0.0 , 0.40150605663 ],
[ -1.21073950301e-06 , 0.0 , 0.40150605663 ],
[ -1.21073950301e-06 , 0.0 , 0.40150605663 ],
[ -1.21073950301e-06 , 0.0 , 0.40150605663 ],
[ -1.21073950301e-06 , 0.0 , 0.40150605663 ],
[ -1.21073950301e-06 , 0.0 , 0.40150605663 ],
[ -1.21073950301e-06 , 0.0 , 0.40150605663 ],
[ -1.21073950301e-06 , 0.0 , 0.40150605663 ],
[ -1.21073950301e-06 , 0.0 , 0.40150605663 ],
[ -1.21073950301e-06 , 0.0 , 0.40150605663 ],
[ -1.21073950301e-06 , 0.0 , 0.40150605663 ],
[ -1.21073950301e-06 , 0.0 , 0.40150605663 ],
[ -1.21073950301e-06 , 0.0 , 0.40150605663 ],
[ -1.21073950301e-06 , 0.0 , 0.40150605663 ],
[ -1.21073950301e-06 , 0.0 , 0.40150605663 ],
[ -1.21073950301e-06 , 0.0 , 0.40150605663 ],
[ -1.21073950301e-06 , 0.0 , 0.40150605663 ],
[ -1.21073950301e-06 , 0.0 , 0.40150605663 ],
[ -1.21073950301e-06 , 0.0 , 0.40150605663 ],
[ -1.21073950301e-06 , 0.0 , 0.40150605663 ],
[ 0.764816906209 , 0.0 , -0.200753505843 ],
[ 0.764816906209 , 0.0 , -0.200753505843 ],
[ 0.764816906209 , 0.0 , -0.200753505843 ],
[ 0.764816906209 , 0.0 , -0.200753505843 ],
[ 0.764816906209 , 0.0 , -0.200753505843 ],
[ 0.764816906209 , 0.0 , -0.200753505843 ],
[ 0.764816906209 , 0.0 , -0.200753505843 ],
[ 0.764816906209 , 0.0 , -0.200753505843 ],
[ 0.764816906209 , 0.0 , -0.200753505843 ],
[ 0.764816906209 , 0.0 , -0.200753505843 ],
[ 0.764816906209 , 0.0 , -0.200753505843 ],
[ 0.764816906209 , 0.0 , -0.200753505843 ],
[ 0.764816906209 , 0.0 , -0.200753505843 ],
[ 0.764816906209 , 0.0 , -0.200753505843 ],
[ 0.764816906209 , 0.0 , -0.200753505843 ],
[ 0.764816906209 , 0.0 , -0.200753505843 ],
[ -0.764815694412 , 0.0 , -0.200752550787 ],
[ -0.764815694412 , 0.0 , -0.200752550787 ],
[ -0.764815694412 , 0.0 , -0.200752550787 ],
[ -0.764815694412 , 0.0 , -0.200752550787 ],
[ -0.764815694412 , 0.0 , -0.200752550787 ],
[ -0.764815694412 , 0.0 , -0.200752550787 ],
[ -0.764815694412 , 0.0 , -0.200752550787 ],
[ -0.764815694412 , 0.0 , -0.200752550787 ],
[ -0.764815694412 , 0.0 , -0.200752550787 ],
[ -0.764815694412 , 0.0 , -0.200752550787 ],
[ -0.764815694412 , 0.0 , -0.200752550787 ],
[ -0.764815694412 , 0.0 , -0.200752550787 ],
[ -0.764815694412 , 0.0 , -0.200752550787 ],
[ -0.764815694412 , 0.0 , -0.200752550787 ],
[ -0.764815694412 , 0.0 , -0.200752550787 ],
[ -0.764815694412 , 0.0 , -0.200752550787 ],
[ -1.21073950301e-06 , 0.0 , 0.40150605663 ],
[ 0.764816906209 , 0.0 , -0.200753505843 ],
[ -0.764815694412 , 0.0 , -0.200752550787 ],
[ -1.21073950301e-06 , 0.0 , 0.40150605663 ],
[ 0.764816906209 , 0.0 , -0.200753505843 ],
[ -0.764815694412 , 0.0 , -0.200752550787 ],
[ 0.764816906209 , 0.0 , -0.200753505843 ],
[ 0.764816906209 , 0.0 , -0.200753505843 ],
[ 0.764816906209 , 0.0 , -0.200753505843 ],
[ 0.764816906209 , 0.0 , -0.200753505843 ],
[ 0.764816906209 , 0.0 , -0.200753505843 ],
[ 0.764816906209 , 0.0 , -0.200753505843 ],
[ 0.764816906209 , 0.0 , -0.200753505843 ],
[ 0.764816906209 , 0.0 , -0.200753505843 ],
[ 0.764816906209 , 0.0 , -0.200753505843 ],
[ -1.21073950301e-06 , 0.0 , 0.40150605663 ],
[ -1.21073950301e-06 , 0.0 , 0.40150605663 ],
[ -1.21073950301e-06 , 0.0 , 0.40150605663 ],
[ -1.21073950301e-06 , 0.0 , 0.40150605663 ],
[ -1.21073950301e-06 , 0.0 , 0.40150605663 ],
[ -1.21073950301e-06 , 0.0 , 0.40150605663 ],
[ -1.21073950301e-06 , 0.0 , 0.40150605663 ],
[ -1.21073950301e-06 , 0.0 , 0.40150605663 ],
[ -0.764815694412 , 0.0 , -0.200752550787 ],
[ -0.764815694412 , 0.0 , -0.200752550787 ],
[ -0.764815694412 , 0.0 , -0.200752550787 ],
[ -0.764815694412 , 0.0 , -0.200752550787 ],
[ -0.764815694412 , 0.0 , -0.200752550787 ],
[ -0.764815694412 , 0.0 , -0.200752550787 ],
[ -0.764815694412 , 0.0 , -0.200752550787 ],
[ -0.764815694412 , 0.0 , -0.200752550787 ],
[ -0.764815694412 , 0.0 , -0.200752550787 ],
[ -1.21073950301e-06 , 0.0 , 0.40150605663 ],
[ -1.21073950301e-06 , 0.0 , 0.40150605663 ],
[ -1.21073950301e-06 , 0.0 , 0.40150605663 ],
[ -1.21073950301e-06 , 0.0 , 0.40150605663 ],
[ -1.21073950301e-06 , 0.0 , 0.40150605663 ],
[ -1.21073950301e-06 , 0.0 , 0.40150605663 ],
[ -1.21073950301e-06 , 0.0 , 0.40150605663 ],
[ -1.21073950301e-06 , 0.0 , 0.40150605663 ],
[ -0.764815694412 , 0.0 , -0.200752550787 ],
[ -0.764815694412 , 0.0 , -0.200752550787 ],
[ 0.764816906209 , 0.0 , -0.200753505843 ],
[ 0.764816906209 , 0.0 , -0.200753505843 ]])
#transpose coordinate array
xyz = np.transpose(coord)
#establish distance matrix
n = len(coord)
dist_vectors = list()
dist_vectors.append(xyz[0])
dist_vectors.append(xyz[1])
dist_vectors.append(xyz[2])
dist_vectors = map(list, zip(*dist_vectors))
d = euclidean_distances(dist_vectors, dist_vectors)
#find two nearest neighbors for each segment
x1 = np.zeros((n))
x2 = np.zeros((n))
y1 = np.zeros((n))
y2 = np.zeros((n))
z1 = np.zeros((n))
z2 = np.zeros((n))
for i in range(n):
x_copy = xyz[0]
y_copy = xyz[1]
z_copy = xyz[2]
d1 = np.delete(d[i], i) #removes distance between segment and itself
x_copy = np.delete(x_copy, i)
y_copy = np.delete(y_copy, i)
z_copy = np.delete(z_copy, i)
j1 = np.argmin(d1) #get indice of minimum distance
x1[i] = x_copy[j1]
y1[i] = y_copy[j1]
z1[i] = z_copy[j1]
d2 = np.delete(d1, j1) #removes minimum distance
x_copy = np.delete(x_copy, j1)
y_copy = np.delete(y_copy, j1)
z_copy = np.delete(z_copy, j1)
j2 = np.argmin(d2) #get indice of second minimum distance
x2[i] = x_copy[j2]
y2[i] = y_copy[j2]
z2[i] = z_copy[j2]
#compute normal vector for each segment based on cross product
normal = list()
forGraphs = list()
for i in range(n):
#make vectors for cross product
v1 = np.zeros((3))
v1[0] = x1[i] - coord[i][0]
v1[1] = y1[i] - coord[i][1]
v1[2] = z1[i] - coord[i][2]
v2 = np.zeros((3))
v2[0] = x2[i] - coord[i][0]
v2[1] = y2[i] - coord[i][1]
v2[2] = z2[i] - coord[i][2]
#make cross product and normalize (normal vector should have a unit norm)
nv = np.cross(v1, v2)
nv = nv / np.linalg.norm(nv)
normal.append(nv)
#check of outwards pointing
atv = np.zeros((3))
atv[0] = atoms[i][0] - coord[i][0]
atv[1] = atoms[i][1] - coord[i][1]
atv[2] = atoms[i][2] - coord[i][2]
th_check = math.acos(np.dot(nv, atv) / (np.linalg.norm(nv) * np.linalg.norm(atv)))
if th_check < (math.pi / 2): #if inwards pointing (i. e. pointing towards underlying atom), normal vector is replaced by its opposite
nv[0] = -nv[0]
nv[1] = -nv[1]
nv[2] = -nv[2]
forGraphs.append(np.array([coord[i][0],coord[i][1],coord[i][2],nv[0],nv[1], nv[2]]))
#plot normal vectors (for checkup)
forGraphs = np.asarray(forGraphs)
X, Y, Z, U, V, W = zip(*forGraphs)
fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
ax.quiver(X, Y, Z, U, V, W)
ax.set_xlim([min(xyz[0])- 1, max(xyz[0]) + 1])
ax.set_ylim([min(xyz[1])- 1, max(xyz[1]) + 1])
ax.set_zlim([min(xyz[2])- 1, max(xyz[2]) + 1])
plt.show()
此代码的前 280 行主要用于重现结果所需的坐标表。这段代码最重要的部分是第 282 行到第 355 行,其中实现了我刚刚概述的算法。
在此先感谢您的帮助!
估计法线的"standard"过程是为每个点找到k个最近的邻居,其中k是一个小数字(十?)。然后计算通过这些点的最佳平面拟合,并使用平面的法线。
不幸的是,法线的符号是不确定的,你需要实施一个法线一致性强制过程。也许在你的情况下这更容易,因为所有法线似乎都指向远离某个中心点。
我按照@Yves Daoust 的建议使用了最近邻技术,并使用奇异值分解 (SVD) 通过 10 个最近邻来拟合一个平面。在视觉上简单的情况下,这仍然会给出错误的答案。我不明白为什么。这是我得到的示例: