网站开发软件开发项目/今日热搜第一名

综述

蓝噪声采样是一种在mesh上重要的采样技术。在迭代次数满足的情况下，他基本满足点的均匀分布。

代码

clc
clear
[V,F]  = readOBJ('baby.obj');
N = random_points_on_mesh(V, F, 4000, 'Color', 'blue', 'MaxIter', 300);
filename = '/Users/frankdura/Desktop/sample.txt';fileID=fopen(filename,'a')fprintf(fileID,'%f %f %f\n',N');fclose(fileID);

核心程序
采样函数
random_points_on_mesh

function [N,I,B,r] = random_points_on_mesh(V,F,n,varargin)% RANDOM_POINTS_ON_MESH Uniform random sampling of a mesh.%% N = random_points_on_mesh(V,F,n)% [N,I,B,r] = random_points_on_mesh(V,F,n,'ParameterName',parameter_value,...)%% Inputs:%   V  #V by dim list of vertex positions%   F  #F by 3 list of triangle indices into V%   n  number of samples%   Optional:%     'Color' followed by "color of noise", one of:%       {'white'}  uniform random noise%       'blue'  approximate blue noise via approximate Poisson disk sampling%     'MaxIter'  followed by maximum number of iterations, only applicable%       for blue noise.%     'Bias'  **experimental** followed by #F list of biases per triangle%       {doublearea(V,F)}.  This bias is a little funny because it effects%       the sampling, but not the distances...probably it gets "washed out"%       by the iterations.% Outputs:%   N  n by dim list of sample positions%   I  n list of face indices upon which N are%   B  n by 3 list of barycentric coordinates so that:%     N =  ...%       bsxfun(@times,B(:,1),V(F(I,1),:)) +  ...%       bsxfun(@times,B(:,2),V(F(I,2),:)) +  ...%       bsxfun(@times,B(:,3),V(F(I,3),:));%   r  radius used for blue noise sampling%% See also: farthest_points%% default valuescolor = 'white';max_iter = 100;A = [];% Map of parameter names to variable namesparams_to_variables = containers.Map( ...{'Color','MaxIter','Bias'}, ...{'color','max_iter','A'});v = 1;while v <= numel(varargin)param_name = varargin{v};if isKey(params_to_variables,param_name)assert(v+1<=numel(varargin));v = v+1;% Trick: use feval on anonymous function to use assignin to this workspace feval(@()assignin('caller',params_to_variables(param_name),varargin{v}));elseerror('Unsupported parameter: %s',varargin{v});endv=v+1;endif isempty(A)A = doublearea(V,F);endassert(size(A,1) == size(F,1));switch colorcase 'white'A = A./sum(A);A0 = [0;A];C = cumsum(A0);R = rand(n,1);[~,I] = histc(R,C);S = rand(n,1);T = rand(n,1);B = [1-sqrt(T) (1-S).*sqrt(T) S.*sqrt(T)];N =  ...bsxfun(@times,B(:,1),V(F(I,1),:)) +  ...bsxfun(@times,B(:,2),V(F(I,2),:)) +  ...bsxfun(@times,B(:,3),V(F(I,3),:));case 'blue'% This is a "cheap hack" way of getting something like Poisson-Disk% sampling which approximates a blue noise sampling.k = 10;%BV = biharmonic_embedding(V,F,k);randsample = @(n) random_points_on_mesh(V,F,n,'Color','white','Bias',A);embed = @(N,I,B) ...bsxfun(@times,B(:,1),V(F(I,1),:)) +  ...bsxfun(@times,B(:,2),V(F(I,2),:)) +  ...bsxfun(@times,B(:,3),V(F(I,3),:));      % initial sampling[N,I,B] = randsample(n);BN = embed(N,I,B);% Expected radius[~,D] = knnsearch(BN,BN,'K',2);D = D(:,2);%r = mean(D);r = mean(D)*sqrt(2);for iter = 1:max_iter% Maybe `rangesearch` would be faster%% Should only query new points...[kI,D] = knnsearch(BN,BN,'K',2);% Discard first match: always selfkI = kI(:,2);D = D(:,2);% Find darts to rethrow, first of each pair with small distancerethrow = unique(kI( kI>(1:numel(kI))' & D<r));%fprintf('iter: %d, #rethrow: %d\n',iter,numel(rethrow));if numel(rethrow) == 0break;end[N(rethrow,:),I(rethrow),B(rethrow,:)] = randsample(numel(rethrow));keep = setdiff(1:n,rethrow);N = N([keep(:);rethrow(:)],:);I = I([keep(:);rethrow(:)],:);B = B([keep(:);rethrow(:)],:);BN = embed(N,I,B);if iter == max_iterwarning('Max iteration (%d) reach without convergence',max_iter);endendendend

readobj

function [V,F,UV,TF,N,NF] = readOBJ(filename,varargin)% READOBJ reads an OBJ file with vertex/face information%% [V,F,UV,TF,N,NF] = readOBJ(filename)% [V,F,UV,TF,N,NF] = readOBJ(filename,'ParameterName',ParameterValue,...)%% Input:%  filename  path to .obj file%  Optional:%    'Quads' whether to output face information in X by 4 matrices (faces%      with degree larger than 4 are still triangulated). A trailing zero%      will mean a triangle was read.% Outputs:%  V  #V by 3 list of vertices%  F  #F by 3 list of triangle indices%  UV  #UV by 2 list of texture coordinates%  TF  #F by 3 list of triangle texture coordinates%  N  #N by 3 list of normals%  NF  #F by 3 list of triangle corner normal indices into N%% Examples:%   % read a quad/triangle mesh and display it%   [V,F] = readOBJ('quads.obj','Quads',true);%   % Turn triangles into degenerate quads %   DF = (F==0).*F(:,[4 1 2 3])+(F~=0).*F;%   trisurf(DF,V(:,1),V(:,2),V(:,3));%%% See also: load_mesh, readOBJfast, readOFF% default valuesquads = false;% Map of parameter names to variable namesparams_to_variables = containers.Map( ...{'Quads'}, ...{'quads'});v = 1;while v <= numel(varargin)param_name = varargin{v};if isKey(params_to_variables,param_name)assert(v+1<=numel(varargin));v = v+1;% Trick: use feval on anonymous function to use assignin to this workspace feval(@()assignin('caller',params_to_variables(param_name),varargin{v}));elseerror('Unsupported parameter: %s',varargin{v});endv=v+1;endnumv = 0;numf = 0;numuv = 0;numtf = 0;numn = 0;numnf = 0;% simplex sizeif quadsss = 4;elsess = 3;end% Amortized array allocationV = zeros(10000,3);F = zeros(10000,ss);UV = zeros(10000,3);TF = zeros(10000,ss);N = zeros(10000,3);NF = zeros(10000,ss);triangulated = false;all_ss = true;fp = fopen( filename, 'r' );type = fscanf( fp, '%s', 1 );count = 0;while strcmp( type, '' ) == 0line = fgets(fp);if strcmp( type, 'v' ) == 1v = sscanf( line, '%lf %lf %lf' );numv = numv+1;if(numv>size(V,1))V = cat(1,V,zeros(10000,3));endV(numv,:) = [v(1:3)'];elseif strcmp( type, 'vt')v = sscanf( line, '%f %f %f' );numuv = numuv+1;if size(UV,2)>2 && length(v) == 2UV = UV(:,1:2);endif(numuv>size(UV,1))UV = cat(1,UV,zeros(10000,length(v)));endUV(numuv,:) = [v'];elseif strcmp( type, 'vn')n = sscanf( line, '%f %f %f' );numn = numn+1;if(numn>size(N,1))N = cat(1,N,zeros(10000,3));endN(numn,:) = [n'];elseif strcmp( type, 'f' ) == 1[t, count] = sscanf(line,'%d/%d/%d %d/%d/%d %d/%d/%d %d/%d/%d %d/%d/%d');if (count>2)tf = t(2:3:end);nf = t(3:3:end);t = t(1:3:end);else[t, count] = sscanf(line, '%d/%d %d/%d %d/%d %d/%d %d/%d');if (count>2)tf = t(2:2:end);t = t(1:2:end);nf = -ones(numel(t),1);else[t, count] = sscanf(line, '%d//%d %d//%d %d//%d %d//%d %d//%d');if (count>2)nf = t(2:2:end);t = t(1:2:end);tf = -ones(numel(t),1);else[t, count] = sscanf( line, '%d %d %d %d %d %d %d %d %d %d %d\n' );if (count>2)tf = -ones(numel(t),1);nf = -ones(numel(t),1);else[t, count] = sscanf( line, '%d// %d// %d// %d// %d// %d// %d// %d// %d// %d// %d//\n' );tf = -ones(numel(t),1);nf = -ones(numel(t),1);endendendendt = t + (t<0).*   (numv+1);tf = tf + (tf<0).*(numuv+1);nf = nf + (nf<0).*(numn+1);assert(numel(t) == numel(tf));assert(numel(t) == numel(nf));if numel(t) > ssif ~triangulatedwarning('Trivially triangulating high degree facets');endtriangulated = true;endj = 2;i = 1;%Vt = V(t,:);%[~,A] = affine_fit(Vt);%VtA = Vt*A;%VtA0 = Vt*A;%[~,alpha] = curvature(VtA);%flip = -sign(sum(alpha));%E = [1:size(VtA,1);2:size(VtA,1) 1]';%[dV,dF] = triangle(VtA,E,[]);%if size(dF,1)>2%  tsurf(dF,dV);%  hold on;%  plot(VtA0([1:end 1],1),VtA0([1:end 1],2),'LineWidth',3);%  hold off%  pause%endwhile trueif numel(t) > sscorners = [1 2 3];%plot(VtA0([1:end 1],1),VtA0([1:end 1],2));%hold on;%plot(VtA([1:3],1),VtA([1:3],2),'LineWidth',3);%hold off;%expand_axis(2);%pause;%[~,alpha] = curvature(VtA,[1 2;2 3]);%alpha = flip * alpha(2);%others = VtA(setdiff(1:end,corners),:);%these = VtA(corners,:);%w = inpolygon(others(:,1),others(:,2),these(:,1),these(:,2));%alpha%if alpha>=0 && ~any(w)%  % lazy%  t = t([2:end 1]);%  VtA = VtA([2:end 1],:);%  continue;%endelseif all_ss && numel(t)<sswarning('Small degree facet found');all_ss = false;endcorners = 1:numel(t);endnumf = numf+1;if(numf>size(F,1))F = cat(1,F,zeros(10000,ss));endF(numf,1:numel(corners)) = [t(corners)'];numtf = numtf+1;if(numtf>size(TF,1))TF = cat(1,TF,zeros(10000,ss));endTF(numtf,1:numel(corners)) = [tf(corners)'];numnf = numnf+1;if(numnf>size(NF,1))NF = cat(1,NF,zeros(10000,ss));endNF(numnf,1:numel(corners)) = [nf(corners)'];if numel(t) <= ssbreak;endt = t([1 3:end]);tf = tf([1 3:end]);nf = nf([1 3:end]);%VtA = VtA([1 3:end],:);if numel(t) < 3break;endendelseif strcmp( type, '#' ) == 1%fscanf( fp, '%s\n', 1 );% ignore lineendtype = fscanf( fp, '%s', 1 );endfclose( fp );%try%    F = cell2mat(F);%endV = V(1:numv,:);F = F(1:numf,:);UV = UV(1:numuv,:);TF = TF(1:numtf,:);N = N(1:numn,:);NF = NF(1:numnf,:);%% transform into array if all faces have the same number of verticesif (size(UV,1)>0) UV = UV; endend