weixin_47616793的博客

原创 matlab提取四面体网格的表面三角网格

function F=extractSurfaceTriangles(T)F_all = listFaces(T);% 提取表面三角网格F_surf = extractSurface(T);F=F_surf(:,1:3);function f = list_faces(t)% t: 四面体网格的面信息，每行代表一个四面体的顶点索引endfunction F_surf = extractSurface(t)f = listFaces(t);f_unique = extractUniqueTrian

原创 用python导入点云和画出一个3D图形的程序

import numpy as npimport open3dimport mayavi.mlabmesh=[[[0]*3]*10000]*14FDI=[17,16,15,14,13,12,11,21,22,23,24,25,26,27];#for i in range(14):#mesh[i]=open3d.io.read_triangle_mesh(‘C:/Users/lixiaohu/Desktop/auto切牙数据/1234567/150451-0/tooth_0_’+str(FDI[i

原创 python读取点云，使用点云，绘3D图，PCA

import numpy as npimport open3dfrom vispy import ioimport numpy as npimport osimport mayavi.mlabmesh=[[[0]*3]*10000]*14FDI=[17,16,15,14,13,12,11,21,22,23,24,25,26,27];for i in range(14):mesh[i]=open3d.io.read_triangle_mesh(‘C:/Users/lixiaohu/Deskt

原创 python怎么批量读取点云，list相当于matlab的cell数组

import numpy as npimport open3dfrom vispy import ioimport numpy as npimport osimport mayavi.mlabmesh=[[[0]*3]*10000]*14FDI=[17,16,15,14,13,12,11,21,22,23,24,25,26,27];for i in range(14):mesh[i]=open3d.io.read_triangle_mesh(‘C:/Users/lixiaohu/Deskt

原创 机器学习的API接口的命令使用方法

这C:\Users\lixiaohu>cd DesktopC:\Users\lixiaohu\Desktop>dir *.zip驱动器 C 中的卷没有标签。卷的序列号是 1A68-AC22C:\Users\lixiaohu\Desktop 的目录2021/10/25 09:25 8,630,059 357176_test.zip2021/05/17 17:07 32,989,448 Netch-完整版.zip2021/07/16 10:32

原创 对icp程序的一些注解，FILEexchange，MATLAB

clear all;syms a b c d e f g h i j k l m n o p q r s t u v w x y z;model=[a b c;d e f;g h i;x y z]data=[k l m n; o p q r; s t u v ]N=4size(data,2)med=mean(data,2);mem=mean(model,1);C=datamodel-(Nmed)mem%这里的矩阵分别为N3，和3*N类型对应于找了最近点之后的两点群的点对%相当于点数不同的

原创 汉明码的解码一例，含四位校验码，matlab

clear alla=[1,0,1,1,0,1,0,1,1,1,0];k=4n=length(a(????); %求输入数组长度h=zeros(1,n-k); %解码和译码的n含义不同check=zeros(1,k); %四个组的验证数组che=0;if k4if a(1)~=rem(a(3)+a(5)+a(7)+a(9)+a(11),2); %四个组，如果有任意一个组验证不通过，则相应数组的变量变为1，标志位che变成1，方便后面修改

原创 汉明码校验位设置后，怎么填入内容呢？

这个图片是一个快速确定归到哪一个类别的汉明码的归类图。汉明码编码时要再源码上插入校验位，在1，2，4，8.，16…等2的幂位置上插入即可。现在问题来了。插了这么多空，怎么填值进去呢？谁对这个值产生贡献呢？归类图就解决这样的问题，例如对9号位置，只需要查表，得到它归为P1和P8类，所以它会为该两位的取值做贡献。怎么填写贡献呢？利用奇偶校验，把归类为该位上的所有数字加起来，根据奇偶性给0或1.就是填写方法。...

原创 汉明码是怎么通过增加校验位来校验的，matlab

cleara=[1,0,1,0,1,1,1];k=4n=length(a(????);z=1;j=1;h=zeros(1,n+k);for i=1:n+k %将输入数组写入全零矩阵待用，且跳过验证if i==2^(z-1) %在2的次幂位置上添加索引，用来填写增加的校验码作为索引的内容z=z+1;continueelseh(i)=a(j) %所以h就是一个膨胀了的数组a，加入了几个校验码j=j+1;endendh =

原创 [7,4]循环码的全部码字用matlab导出，并验证生成多项式的理想结论

n=7;k=4;p=cyclpoly(n,k,‘all’);[H,G]=cyclgen(n,p(1,:));a=G(1,:)b=G(2,:)c=G(3,:)d=G(4,:)a+ba+ca+db+cb+dc+da+b+cb+c+da+c+da+b+da+b+c+d运行如下： 1 0 1 1 0 0 0 1 1 1 0 1 0 0 1 1 0 0

原创 preall点云粗对齐作为matlab函数，pca和普氏分析法都有函数，主成分也有函数

target=V0;十万点source=V111;七万点[COEFF,Prealligned_source] = princomp(source);[COEFF,Prealligned_target] = princomp(target);% the direction of the axes is than evaluated and corrected if necesarry.Maxtarget=max(Prealligned_source)-min(Prealligned_source);

原创 极简的有限差分法做ICP的主干matlab代码

极简的有限差分法做ICP的主干代码npoinst=10000;x=rand(npoinst,1)*100-50; y=rand(npoinst,1)*100-50; z=sqrt(x.2+y.2);Points_Static=[x y z]; x=rand(npoinst-100,1)*100-50; y=rand(npoinst-100,1)*100-50; z=sqrt(x.2+y.2);Points_Moving=[x y z];scale_matrix = zeros(4,4);defa

原创 普氏分析法的计算过程，使用的误差，matlab

syms a b c d e f A B C D E FX=[a,b;c,d;e,f];Y=[A,B;C,D;E,F];% X=[1,2;3,3;4,5]% Y=[1,2;3,4;4,2]muX = mean(X,1)muY = mean(Y,1)X0 = X - repmat(muX, 3, 1)Y0 = Y - repmat(muY, 3, 1)ssqX = sum(X0.^2,1)ssqY = sum(Y0.^2,1)ssqX = sum(ssqX)ssqY = sum(ssq

原创 普氏分析法的输入和输出实例，matlab图像处理

X=[1,2;3,3]Y=[1,2;3,4][a,b,c]=procrustes(X,Y)Z = c.b * Y * c.T + c.ca =-4.4409e-16b = 1 2 3 3c =T: [2x2 double]b: 0.7906c: [2x2 double]Z =1.0000 2.00003.0000 3.0000c.Tc.bc.cans =旋转矩阵0.3162 0.94870.9487 -0.3162a

原创 两个点集的距离，SVD的实施方法

A = X0’ * Y0;[L, D, M] = svd(A);T = M * L’;if isempty(doReflection) % ‘best’% Let the data decide if a reflection is needed.elsehaveReflection = (det(T) < 0);% If we don’t have what was asked for …if (doReflection ~= haveReflection)% … then eit

原创 normrow行范数，matlab

function [ B ] = normrow( A )% NORMROW Compute l2 row vector norms%% B = normrow( A )%% Input:% A #A by D list of row vectors of dimension D% Output:% B #A list of norms of row vectors in A%% Copyright 2011, Alec Jacobson (jacobson@inf.ethz.

原创 点云配准思路

V0是原始点云，V1是修剪后的点云，V0经过旋转得到V111，快要与V1才走到中途V111处就不动了，接下来精细化调整向V0。我们从V111上取点P，计算它离V0的哪一个点最近。发现是第79799个点离它最近。接下来。取V0（79799）得点PP...

原创 random的用法，matlab，正态泊松均匀分布

泊松分布K=[1,2,3,4]random(‘poisson’,K )得到K = 1 2 3 4ans = 2 0 2 2取了四个数，为什么是四个？K=[1,2,3,4]random(‘poisson’,1,K )ans(:,:,1,1) = 2 1ans(:,:,2,1) = 2 1ans(:,:,3,1) = 0 1ans(:,:,1,2) = 1 2ans(:,:,2,2)

原创 找到第一个不平凡的整基的三次扩张$x^3+3x+5$

def f(p,q) :K = NumberField(x^3 +p*x+q, ‘a’);return K.integral_basis()[print(p,q,f(p,q)) for p in range(1,7,2) for q in range(1,7,2)]1 1 [1, a, a^2]1 3 [1, a, a^2]1 5 [1, a, a^2]3 1 [1, a, a^2]3 3 [1, a, a^2]35[1/3∗a2+1/3∗a+1/3,a,a2]3 5 [1/3*a^2 +

原创 sagemath写的实数展成Hirzebruch连分数与高斯二次域的类数

def f®:d=[0]*200e=[0]*200d[0]=ceil®;e[0]=1/(d[0]-r);for i in range(1,8):d[i]=ceil(e[i-1]);e[i]=1/(d[i]-e[i-1]);M=MatrixSpace(ZZ,8,1)d=M(d)return transpose(d)f(N(sqrt(163)))二次域得到类数与上面的连分数有关，周期为n1，。。。，nr时，其和/3-r为二次域的类数，在p为模4余3时成立def f§:K.=Quadr

原创 费马三次曲面的有理参数化

x3+y3+z3+w3=0=加减换元后为:X3+3XY2+Z3+3ZW2=0x^3+y ^3+z ^3 +w ^3=0\\\\\\=加减换元后为:\\X ^3+3X Y ^2+Z ^3+3ZW ^2=0x3+y3+z3+w3=0=加减换元后为:X3+3XY2+Z3+3ZW2=0仿射化后:1+3Y2+Z3+3ZW2=0仿射化后:\\1+3Y ^2+Z ^3 +3Z W ^2=0仿射化后:1+3Y2+Z3+3ZW2=0这是一个带参数W的椭圆曲线。怎么找出他的有理参数化？sagemath代码：va

原创 laplace变形系统重建时，固定控制点在变换后，移动了，mathematica

v = {{Subscript[x, 1], Subscript[y, 1]}, {Subscript[x, 2], Subscript[y, 2]}, {Subscript[x, 3], Subscript[y, 3]}, {Subscript[x, 4],Subscript[y, 4]}, {Subscript[x, 5], Subscript[y, 5]}, {Subscript[x,6], Subscript[y, 6]}}A = {{1, -1/3, 0, -1/3, -1/3, 0},

原创 laplace形变用于数学系统构建，图像变形算法，Matlab

clear allload seamountN = 6;v=[1 0.5;0.78023 0.44678;0.081126 0.30635;0.6 0.2;0.48679 0.81763;0.2586 0.61483]tri= delaunay(v(:,1),v(:,2))z=zeros(N,1);% Delaunay triangulationfigure;subplot(2,1,1);triplot(tri,v(:,1

原创 拉普拉斯坐标（微分坐标）怎么求，laplace算法图像变形，matlab

clear allload seamountN = 8;x = rand(N,1);y = rand(N,1);z=zeros(N,1);% Delaunay triangulationtri = delaunay(x,y)DT = delaunayTriangulation(x,y)% Calculate adjacency matrixAdjMat = false(N);for kk = 1:size(tri,1)AdjMat(tri(kk,1), tri(kk,2)) = tr

原创 Delaunay triangulation与delaunay有什么不同，前者是类，后者是矩阵，matlab

load seamountN = 8;x = rand(N,1);y = rand(N,1);z=zeros(N,1);% Delaunay triangulationtri = delaunay(x,y)DT = delaunayTriangulation(x,y)tri = 1 3 8 1 4 3 1 2 6 2 5 7 6 2 7 1 8 2 3 4 8 8

原创 matlab用适用于三维的trisurf画二维三角剖分，让z=0

load seamountN = 8;x = rand(N,1);y = rand(N,1);z=zeros(N,1);% Delaunay triangulationtri = delaunay(x,y);% DT = delaunayTriangulation(x,y)trisurf(tri,x,y,z)

原创 matlab二维点集的德洛内三角化的邻接矩阵，triplot画三角网格

N = 8;x = rand(N,1);y = rand(N,1);% Delaunay triangulationtri = delaunay(x,y);DT = delaunayTriangulation(x,y)% Calculate adjacency matrixAdjMat = false(N);for kk = 1:size(tri,1)AdjMat(tri(kk,1), tri(kk,2)) = true;AdjMat(tri(kk,2), tri(kk,3)) = tr

原创 matlab程序的函数太多的管理技巧

matlab程序的函数太多了，一般点检一下之后，每次上班都像点检消火栓一样点检一下之后，可以避免无所事事的发呆。一般放了一个长假过后就是这种感觉。每个matlab的函数是干嘛都不知道了。...

原创 mathematica怎么查找和筛选矩阵的元素

f[k_] := k[[1]] + k[[2]]Select[{{1, 2}, {4, 7}, {6, 2}}, f[#] > 5 &]运行结果：{{4, 7}, {6, 2}}

原创 Laplace图像变形的余切矩阵的效果，Matlab

clear all;mesh_source = ‘C:\Users\li\Desktop\101028\101028_l\T_0.obj’;[v,F] = load_mesh(mesh_source);cotmatrix(v,F)结果：非对角矩阵。对角元+其他零散的元素(5158,5227) 0.20947(5159,5227) 1.06(5160,5227) 0.45512(5226,5227) 0.93866(5227,5227)

原创 Mathematica求出mass matrix用于双调和图像变形，类型为barycentric型的

ClearAll[“Global`*”]a = Import[“T.obj”, “Table”]f1 = Coefficient[First[Total[a]], “f”, 1]a1 = Coefficient[First[Total[a]], “v”, 1]V = a[[1 ;; a1, 2 ;; 4]]F = a[[a1 + 1 ;; f1 + a1, 2 ;; 4]]F = Transpose[F];i1 = F[[1]];i2 = F[[2]];i3 = F[[3]];v1 =

原创 matlab的size和sum对行列求和的视觉设计记忆

原创 mass矩阵入门，matlab

clear all;mesh_source = ‘C:\Users\lixiaohu\Desktop\101028\101028_l\TOOTH_0.obj’;[v,F] = load_mesh(mesh_source);massmatrix(v,F,‘barycentric’)(4592,4592) 0.033451(4593,4593) 0.018297(4594,4594) 0.059208(4595,4595) 0.126(4596,4596)

原创 matlab三维对称矩阵高速计算求逆矩阵，先用符号验证一遍。

syms a b c d e f;Ia = (d.*f-e.^2)Ib = (c.*e-b.*f)Id =(a.*f-c.^2)Ic= (b.*e-c.*d)Ie=(b.*c-a.*e)If=(a.*d-b.^2)[Ia,Ib,Ic;Ib,Id,Ie;Ic,Ie,If]expand([Ia,Ib,Ic;Ib,Id,Ie;Ic,Ie,If][a ,b,c;b,d,e;c,e f])/ (- fb^2 + 2bce - dc^2 - ae^2 + ad*f)

原创 移动最小二乘MLSD2DTransform函数的结构，Matlab

原创 mathmatica矩阵的运算，相乘，转置，求逆矩阵

原创 最小二乘法时，对三维的点集会遇到解冗余的方程AX=b的问题。 此时A‘AX=A‘b的解法很实用。但是A’A的逆矩阵，maple

计算最小二乘法时，对三维的点集用直线拟合。会遇到解冗余的方程AX=b的问题。此时A‘AX=A’b的解法很实用。但是A’A的逆矩阵，编程难。备忘其逻辑结构如下。

原创 移动最小二乘法MLSD2D的数据预先处理的函数PrecomputeAffine的数学逻辑，matlab

function data = PrecomputeAffine(p,v,w)Pstar = PrecomputeWCentroids(p,w);M1 = v - Pstar;a = zeros(1,size(Pstar,2));b = a;d = a;Phat = cell(1,size(p,2));for i=1:size(p,2)Phat{i} = repmat(p(:,i),[1,size(Pstar,2)])-Pstar; a = a + w(i,:).*Phat{i}(1

原创 图像变形MLSD2DTransform---matlab自带的函数，其中内嵌的子函数PrecomputeAffine预先计算了权重，对三维数据不可！要改造

clear allsyms a1 a2 a3 a4 b1 b2 b3 b4 d1 d2 d3 d4 M1;a=[a1,a2,a3,a4];b=[b1,b2,b3,b4];d=[d1,d2,d3,d4];det = a.*d - b.^2//det会是一个向量，每个元都是a1d1-b1b1这样类似的元。Ia = d./detIb = -b./detId = a./detF1 = [sum(M1.[Ia;Ib],1);sum(M1.[Ib;Id],1)]运行结果---------------

原创 图像变形MLSD2DTransform---matlab自带的函数，其中内嵌的子函数PrecomputeWeights预先计算了权重，对三维数据可以使用

clear all;mesh_source = ‘C:\Users\lixiaohu\Desktop\101028\101028_l\TOOTH_0.obj’;[v,F] = load_mesh(mesh_source);v=v(1:16,:)v=v’size(v,1)p = [-1.0843 -1.2803 -1.5304 -1.4777 ;-3.1908 -3.1335 -3.0674 -3.0253;-3.1908