From 90b2d9ba092d3cd7c4d621ade8667182593af876 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?=E5=88=98=E6=98=8E=E5=AE=8F?= Date: Thu, 22 Jun 2023 06:27:19 +0000 Subject: [PATCH] Code library MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Signed-off-by: 刘明宏 --- lib/AddBlocks.m | 33 +++ lib/AddRecParametricSurface.m | 34 +++ lib/AndInterpolationCurve.m | 13 + lib/AndPoint.m | 13 + lib/Comsol_with_Matlab_Start.m | 19 ++ lib/Variables 1.txt | 9 + lib/Variables 2.txt | 12 + lib/get_curveTxtFile.m | 17 ++ lib/interpclosed.m | 420 +++++++++++++++++++++++++++++++++ lib/lofting.m | 30 +++ lib/plotSlice.m | 27 +++ lib/smooth1q.m | 241 +++++++++++++++++++ lib/smooth1qExample.m | 30 +++ 13 files changed, 898 insertions(+) create mode 100644 lib/AddBlocks.m create mode 100644 lib/AddRecParametricSurface.m create mode 100644 lib/AndInterpolationCurve.m create mode 100644 lib/AndPoint.m create mode 100644 lib/Comsol_with_Matlab_Start.m create mode 100644 lib/Variables 1.txt create mode 100644 lib/Variables 2.txt create mode 100644 lib/get_curveTxtFile.m create mode 100644 lib/interpclosed.m create mode 100644 lib/lofting.m create mode 100644 lib/plotSlice.m create mode 100644 lib/smooth1q.m create mode 100644 lib/smooth1qExample.m diff --git a/lib/AddBlocks.m b/lib/AddBlocks.m new file mode 100644 index 0000000..d8a6acb --- /dev/null +++ b/lib/AddBlocks.m @@ -0,0 +1,33 @@ +function [object_blk,blk_pos_corner ]=AddBlocks(model,blk) +%Add blocks, used in model expansion and loading rectangular sources on the topographic surface + +blk_pos_corner = []; +ins0 = 0; +for n = 1:size(blk.Lx,2) + %Size unit：m + lenx = blk.Lx(n); + leny = blk.Ly(n); + lenz = blk.Lz(n); + %Center position coordinate + xc = blk.CentCoord(n,1); + yc = blk.CentCoord(n,2); + zc = blk.CentCoord(n,3); + + blk_position = [xc yc zc]; + blk_size = [lenx leny lenz]; + blkLabel = ['blk' num2str(n+ins0)]; + + model.geom('geom1').feature.create(blkLabel,'Block'); + model.geom('geom1').feature(blkLabel).set('size',blk_size); + model.geom('geom1').feature(blkLabel).set('pos',blk_position); + model.geom("geom1").feature(blkLabel).set("rot", blk.angle); + model.component("mod1").geom("geom1").feature( blkLabel).set("base", "center"); + + object_blk{n} = blkLabel; + temp = []; + blk_pos_corner = cat(1, blk_pos_corner, temp); +end + +model.component("mod1").geom("geom1").run(); + + diff --git a/lib/AddRecParametricSurface.m b/lib/AddRecParametricSurface.m new file mode 100644 index 0000000..5b546dc --- /dev/null +++ b/lib/AddRecParametricSurface.m @@ -0,0 +1,34 @@ +function [object_rec,ps_pos_rec ]=AddRecParametricSurface(model,rec,lengthcurve,heightcurve) +% +%A vertical auxiliary parametric surfaces are added for constructing the +% receives on the terrain surface +% ins0 = 0; +% xrec=[2 10]; +% yrec=[0 0]; +% zrec=[0 0]; +% lengthcurve=[5 5]; +% heightcurve=5; + +x = rec(:,1); +y = rec(:,2); +z = rec(:,3); +nrec = length(x); +object_rec=cell(1,nrec);%{}; +ps_pos_rec=[]; +ins0 = 0; +for n = 1:nrec + psLabel = ['ps' num2str(n+ins0)]; + ps(n+ins0)= model.geom('geom1').feature.create( psLabel ,'ParametricSurface'); + model.geom('geom1').feature( psLabel ).set('parmin1',num2str(x(n)-lengthcurve/2)); + model.geom('geom1').feature( psLabel ).set('parmin2',num2str(z(n)-heightcurve)); + model.geom('geom1').feature( psLabel ).set('parmax1',num2str(x(n)+lengthcurve/2)); + model.geom('geom1').feature( psLabel ).set('parmax2',num2str(z(n)+heightcurve)); + model.geom('geom1').feature( psLabel ).set('coord',{'s1',num2str(y(n)),'s2'}); + model.geom('geom1').feature( psLabel ).set('maxknots',{'4'}); + object_rec{n} = psLabel; + ps_pos_temp = [x(n)-lengthcurve/2 y(n) z(n)-heightcurve x(n)-lengthcurve/2 y(n) z(n)+heightcurve]; + ps_pos_rec = cat(1, ps_pos_rec, ps_pos_temp); +end +model.component("mod1").geom("geom1").run(); + + diff --git a/lib/AndInterpolationCurve.m b/lib/AndInterpolationCurve.m new file mode 100644 index 0000000..4f601cd --- /dev/null +++ b/lib/AndInterpolationCurve.m @@ -0,0 +1,13 @@ +function objectIPC = AndInterpolationCurve(model,filenames) +%Interpolating curves are added, which can be used to construct 3D irregular volumes +for n = 1:size(filenames,1) + IPCname =['ipc' num2str(n)]; + model.geom('geom1').create(IPCname, 'InterpolationCurve'); + model.geom('geom1').feature(IPCname).set('type', 'closed'); + model.geom('geom1').feature(IPCname).set('source', 'file'); + model.geom('geom1').feature(IPCname).set('filename', filenames{n}); + model.geom('geom1').feature(IPCname).set('struct', 'sectionwise'); + objectIPC{n} = IPCname; +end +model.geom('geom1').run; + diff --git a/lib/AndPoint.m b/lib/AndPoint.m new file mode 100644 index 0000000..f7fb786 --- /dev/null +++ b/lib/AndPoint.m @@ -0,0 +1,13 @@ + +function objectIPC = AndPoint(model,p) +% Add points that can be used for the endpoints of 3D irregular volumes +for np = 1:size(p,1) + PTname =['pt' num2str(n)]; + model.geom("geom1").create(PTname, "Point"); + model.geom("geom1").feature(PTname).setIndex("p", p(np,1), 0); + model.geom("geom1").feature(PTname).setIndex("p", p(np,2), 1); + model.geom("geom1").feature(PTname).setIndex("p", p(np,3), 2); + objectPT{n} = PTname; +end +model.geom('geom1').run; +end diff --git a/lib/Comsol_with_Matlab_Start.m b/lib/Comsol_with_Matlab_Start.m new file mode 100644 index 0000000..0b59ab6 --- /dev/null +++ b/lib/Comsol_with_Matlab_Start.m @@ -0,0 +1,19 @@ + + +% Start the 'COMSOL Multiphysics with MATLAB' interfaces +% Required: '...\comsolmphserver.exe' and '...\Multiphysics\mli' file address +% You can also manually launch the executable COMSOL Multiphysics with MATLAB. exe + +path = pwd; + +try + mphtags -show + warning('Already connected to a server!'); +catch + winopen('D:\Software\COMSOL60\Multiphysics\bin\win64\comsolmphserver.exe'); + cd 'D:\Software\COMSOL60\Multiphysics\mli'; + mphstart; +end +cd(path); + + diff --git a/lib/Variables 1.txt b/lib/Variables 1.txt new file mode 100644 index 0000000..8e7d0a8 --- /dev/null +++ b/lib/Variables 1.txt @@ -0,0 +1,9 @@ +rho_xy ((abs((Ex2_G*mf.Hx-Ex1_G*mf2.Hx)/(mf.Hx*mf2.Hy-mf2.Hx*mf.Hy)))^2/(2*pi*freq*mu0_const)) "Apparent resistivity, xy" +rho_yx ((abs((Ey1_G*mf2.Hy-Ey2_G*mf.Hy)/(mf.Hx*mf2.Hy-mf2.Hx*mf.Hy)))^2/(2*pi*freq*mu0_const)) "Apparent resistivity, yx" +rho_xx ((abs((Ex1_G*mf2.Hy-Ex2_G*mf.Hy)/(mf.Hx*mf2.Hy-mf2.Hx*mf.Hy)))^2/(2*pi*freq*mu0_const)) "Apparent resistivity, xx" +rho_yy ((abs((Ey2_G*mf.Hx-Ey1_G*mf2.Hx)/(mf.Hx*mf2.Hy-mf2.Hx*mf.Hy)))^2/(2*pi*freq*mu0_const)) "Apparent resistivity, yy" +phi_xy arg(1[S]*(Ex2_G*mf.Hx-Ex1_G*mf2.Hx)/(mf.Hx*mf2.Hy-mf2.Hx*mf.Hy))[rad] "Apparent resistivity phase, xy" +phi_yx arg(1[S]*(Ey1_G*mf2.Hy-Ey2_G*mf.Hy)/(mf.Hx*mf2.Hy-mf2.Hx*mf.Hy))[rad] "Apparent resistivity phase, yx" +phi_xx arg(1[S]*(Ex1_G*mf2.Hy-Ex2_G*mf.Hy)/(mf.Hx*mf2.Hy-mf2.Hx*mf.Hy))[rad] "Apparent resistivity phase, xx" +phi_yy arg(1[S]*(Ey2_G*mf.Hx-Ey1_G*mf2.Hx)/(mf.Hx*mf2.Hy-mf2.Hx*mf.Hy))[rad] "Apparent resistivity phase, yy" +S abs((mf2.Ex/mf2.Hx+mf.Ey/mf.Hy)/(mf.Ex/mf.Hy-mf2.Ey/mf2.Hx)) "" diff --git a/lib/Variables 2.txt b/lib/Variables 2.txt new file mode 100644 index 0000000..9a2f371 --- /dev/null +++ b/lib/Variables 2.txt @@ -0,0 +1,12 @@ +Ex1_G real(mf.Ex)+mf.omega*1[s]*imag(d(mf.psi0,x))/1[S/m]+i*(imag(mf.Ex)+mf.omega*1[s]*real(d(mf.psi0,x))/1[S/m]) "" +Ey1_G real(mf.Ey)+mf.omega*1[s]*imag(d(mf.psi0,y))/1[S/m]+i*(imag(mf.Ey)+mf.omega*1[s]*real(d(mf.psi0,y))/1[S/m]) "" +Ez1_G real(mf.Ez)+mf.omega*1[s]*imag(d(mf.psi0,z))/1[S/m]+i*(imag(mf.Ez)+mf.omega*1[s]*real(d(mf.psi0,z))/1[S/m]) "" +Ex2_G real(mf2.Ex)+mf2.omega*1[s]*imag(d(mf2.psi0,x))/1[S/m]+i*(imag(mf2.Ex)+mf2.omega*1[s]*real(d(mf2.psi0,x))/1[S/m]) "" +Ey2_G real(mf2.Ey)+mf2.omega*1[s]*imag(d(mf2.psi0,y))/1[S/m]+i*(imag(mf2.Ey)+mf2.omega*1[s]*real(d(mf2.psi0,y))/1[S/m]) "" +Ez2_G real(mf2.Ez)+mf2.omega*1[s]*imag(d(mf2.psi0,z))/1[S/m]+i*(imag(mf2.Ez)+mf2.omega*1[s]*real(d(mf2.psi0,z))/1[S/m]) "" +normE1_G sqrt(Ex1_G^2+Ey1_G^2+Ez1_G^2) "" +normE2_G sqrt(Ex2_G^2+Ey2_G^2+Ez2_G^2) "" +normE2_G2 sqrt(Ey2_G^2+Ez2_G^2) "" +Eyz_r (sqrt(real(Ey2_G)^2+real(Ez2_G)^2)) "" +Eyz_i (sqrt(imag(Ey2_G)^2+imag(Ez2_G)^2)) "" +Exy_i (sqrt(imag(Ex2_G)^2+imag(Ey2_G)^2)) "" diff --git a/lib/get_curveTxtFile.m b/lib/get_curveTxtFile.m new file mode 100644 index 0000000..da0acad --- /dev/null +++ b/lib/get_curveTxtFile.m @@ -0,0 +1,17 @@ + +%Gets the path of a file with a specific character in the folder +function filenames = get_curveTxtFile(data_dir,id,str) +% data_dir：folder path +% id：The position of the character +% Finding characters +% example: curve_01.txt;curve_02.txt +% CurveFiles = get_curveTxtFile(data_dir,[1:5],'curve'); + +D = dir(data_dir ); +nf = 0; +for i=3:numel(D) + if strcmp(D(i).name(id),str) + nf = nf+1; + filenames{nf,1} = [data_dir '\' D(i).name] ; + end +end \ No newline at end of file diff --git a/lib/interpclosed.m b/lib/interpclosed.m new file mode 100644 index 0000000..6e5d6a1 --- /dev/null +++ b/lib/interpclosed.m @@ -0,0 +1,420 @@ +function varargout = interpclosed(x,y,varargin) +% INTERPCLOSED Arc-length interpolation, perimeter and area of 2D closed curves defined by points +% +% xyq = INTERPCLOSED(x,y,tq) Interpolates new data points xyq at given +% positions defined by an arc-length parametrization and the query points +% tq, along the closed curve defined by the points specified by x and y. +% The default method 'spline' is used. tq is a n-by-1 array with its +% elements constrained within [0,1], with 0 being the first point of the +% curve and 1 the last one. +% +% [len,area] = INTERPCLOSED(x,y) If tq is not specified and two output +% variables are requested, then only the perimeter and area of the +% interpolated curve are returned. Both outputs are obviously scalars. +% +% [len,area,c] = INTERPCLOSED(x,y) If tq is not specified and three output +% variables are requested, then the perimeter, area and centroid of the +% interpolated curve are returned. The centroid is a vector with the +% position as (x,y). +% +% [len,area,c,I] = INTERPCLOSED(x,y) If tq is not specified and four output +% variables are requested, then the perimeter, area, centroid, and second +% moments of area of the interpolated curve are returned. The second moment +% of area is a vector with three elements: (Ixx, Iyy, Ixy). +% +% pp = INTERPCLOSED(x,y,'pp') If only one output argument is defined and +% the string 'pp' is given as input, the returned variable is the piesewise +% polynomial pp, regardless of the definition of tq. +% +% [xyq,len,area] = INTERPCLOSED(x,y,tq) If tq is provided and there are +% three output variables, the perimeter and area are additionally returned. +% +% [xyq,len,area,c,I] = INTERPCLOSED(x,y,tq) If tq is provided and there are +% five output variables, the perimeter, the area, the centroid, and the +% second moments of area are additionally returned. +% +% [___] = INTERPCLOSED(___,method) By specifying the string method it +% is possible to change how the interpolated points are obtained. +% 'linear': Linear interpolation. The interpolated value at a query point +% is based on linear interpolation of the values at neighboring +% points in each respective dimension. This is the fastest +% method. +% 'spline': Spline interpolation using periodic end conditions. The +% interpolated value at a query point is based on a cubic +% interpolation of the values at neighboring points in each +% respective dimension. This is the default method. +% 'pchip': Shape-preserving piecewise cubic interpolation. The +% interpolated value at a query point is based on a shape- +% preserving piecewise cubic periodic interpolation of the +% values at neighboring points. +% +% [___] = INTERPCLOSED(___,print) By setting the boolean variable print to +% true, more details about the interpolation can be obtained. +% +% Examples: +% 1) Use the spline method to interpolate at 64 query points within the +% 8 points used to sample the original circle: +% +% t = linspace(0,2*pi,9); +% x = sin(t) + 0.2; y = cos(t) - 0.2; +% [len,area,c,I] = interpclosed(x,y); +% fprintf(['Perimeter: %4.5f, Area: %4.5f,\n',... +% 'Centroid: [%4.5f %4.5f], Iz: %4.5f\n',... +% 'To be compared to\n',... +% '2*pi: %4.5f, pi: %4.5f,\n',... +% 'Centroid: [%4.5f %4.5f], pi/2: %4.5f\n'],... +% len,area,c,(I(1)-area*c(2)^2+I(2)-area*c(1)^2),... +% 2*pi,pi,[0.2,-0.2],pi/2) +% +% 2) Get the piecewise polyonomial of a linear interpolation of a given +% set of points and use the polynomial to make a plot: +% +% x = [0 .82 .92 0 -.92 -.82]; y = [.66 .9 0 -.83 0 .9]; +% pp = interpclosed(x,y,'pp','linear'); +% tq = min(pp.breaks):0.001:max(pp.breaks); +% xyq = ppval(pp,tq); +% figure, plot(xyq(1,:),xyq(2,:)) +% +% Find more examples in the File Exchange website. +% +% See also CSCVN, PCHIP, MKPP, PPVAL, INTERPARC, ARCLENGTH, INTERP1. +% +% Author: Santiago M. Benito +% Ruhr-Universität Bochum +% ------------------------------------------------------------------------- +% Contact: santiago.benito@rub.de +% ------------------------------------------------------------------------- +% Current version: 3.0 +% ------------------------------------------------------------------------- +% Last updated: 17.05.2021 +% Changes: +% - It is now possible to compute the second moments of area of the +% fit with this function. + +%% Manage input, output and catch eventual problems. +% Check for errors in the given inputs. +if nargin < 2 + error('INTERPCLOSED:insufficientarguments', ... + 'At least x and y must be supplied.') +end + +if ~isvector(x) || ~isvector(y) || (length(x) ~= length(y)) + error('INTERPCLOSED:baddimension', ... + 'x and y must be vectors of the same length.') +end + +% Set defaults. +method = 'spline'; +print = false; +geomcalc = false; tqgiven = false; +pp = false; + +% Initialize output variables. +len = 0; +area = 0; +c = zeros(1,2); +Ixx = 0; Iyy = 0; Ixy = 0; + +% Check for other input arguments. +if numel(varargin) > 0 + % At least one other argument was supplied. + for ii = 1:numel(varargin) + arg = varargin{ii}; + if ischar(arg) + % It can be the method or the 'pp'-flag. + validstrings = {'pp','linear' 'pchip' 'spline'}; + ind = strncmp(arg,validstrings,2); + if isempty(ind) || (sum(ind) == 0) || (sum(ind) > 1) + error('INTERPCLOSED:invalidmethod', ... + ['Invalid method indicated. Only ''linear'',',... + '''pchip'',''spline'' allowed.']) + end + if ind(1) == 1 + pp = true; + else + method = validstrings{ind>0}; + end + elseif islogical(arg) + % It must be the print variable, set the print sampling distance. + if ~tqgiven, tq = 0:1/32:1; end + print = arg; + else + % It must be tq, defining the parametric arc-length query + % points + if ~isnumeric(arg) + error('INTERPCLOSED:badtq', ... + 'tq must be numeric.') + else + if max(arg) > 1 || min(arg) < 0 + error('INTERPCLOSED:badtq', ... + 'tq elements must be bigger than 0 and smaller than 1.') + end + tqgiven = true; + tq = arg; + end + end + end +end + +% If only one output variable is requested and the 'pp' flag was given, no +% need to compute the interpolations, regardless of the definition of tq. +% If three are given, geometry computations will be needed. If two or three +% are given, but tq was not provided, also compute the geometry computations. +if nargout == 1 && pp && ~print + tqgiven = false; +elseif (nargout == 2 || nargout == 3 || nargout == 4) && ~tqgiven + geomcalc = true; +elseif (nargout == 3 || nargout == 4 || nargout == 5) + geomcalc = true; + if ~tqgiven + error('INTERPCLOSED:badtq', ... + 'tq was not defined and is needed for interpolation.') + end +end + +% Be sure everything is formatted correclty and group it. +x = x(:)'; y = y(:)'; +points = [x;y]; + +% Round to the 15th decimal position to avoid rounding errors. This is +% necessary for the function to recongnize start and ending points +% properly. +points = round(points,15); + +% If the set of points does not describe a closed loop, close it. +if sum(points(:,1) ~= points(:,end)) > 0 + points = [points,points(:,1)]; +end + +% If less than three distinct points are given, no closed curve can be +% formed. +d = sum((diff(points.').^2).'); +if numel(x) - sum(d==0)-1 < 2 + error('INTERPCLOSED:baddimension', ... + 'x and y must be vectors describing at least three distinct points.') +end + +%% Actual program start +% Compute the coefficients of the fit-polynomials according to the user's +% choice. +if strcmpi(method,'linear') + % Remove segments with length equal to zero, the linear interpolation + % has no continuous derivatives anyway. + points(:,d==0) = []; + + % Compute the linear coefficients of the parametric versions of the + % lines. First compute the lengths of each segment, then the cumulative + % length and finally use the slope in each direction to get the coefs. + seglen = sqrt(sum(diff(points,[],2).^2,1)); + cumarc = [0,cumsum(seglen)]; + coefX = [diff(points(1,:))./diff(cumarc);points(1,1:(end-1))]; + coefY = [diff(points(2,:))./diff(cumarc);points(2,1:(end-1))]; + + % Create a piecewise polynomial with the given coefficients. + coefs = zeros(size([coefX,coefY])); + coefs(:,1:2:end) = coefX; + coefs(:,2:2:end) = coefY; + curve = mkpp(cumarc,coefs',2); + + % Provide the differentiation array for later use. + diffarray = [0 0 1;0 0 0]; + + % Since we already have the lenghts of the individual segments, just + % sum everything up and save some time. + len = sum(seglen); + +elseif strcmpi(method,'spline') + % MATLAB(R) already has a very useful function that makes all the work + % for us. + curve = cscvn(points); + + % Provide the differentiation array for later use. + diffarray = [3 0 0;0 2 0;0 0 1;0 0 0]; + +elseif strcmpi(method,'pchip') + % Like in the function CSCVN, if the user specified a point where the + % 2nd derivative is equal to zero, we have to be able to handle the + % situation. + d = sum((diff(points.').^2).'); + + if all(d > 0) + % The fit is periodic. To have the start and end slopes equal to + % each other, some tricks must be done. Extra points will be added + % right before the start and right after the end of the data set. + % The fit will be performed with these points, and then the extra + % pieces will be removed from the general fit. + %pointsNew = [x(end-2:end-1),x,x(2:3);y(end-2:end-1),y,y(2:3)]; + pointsNew = [points(:,end-2:end-1),points,points(:,2:3)]; + + % We need the arc length of the modified dataset, therefore we will + % compute it here. + seglen = sqrt(sum(diff(pointsNew,[],2).^2,1)); + cumarc = [0,cumsum(seglen)]; + + % Fit coefficients are obtained from the MATLAB(R) original pchip + % function. + temp = pchip(cumarc,pointsNew(1,:)); coefX = temp.coefs; + temp = pchip(cumarc,pointsNew(2,:)); coefY = temp.coefs; + + % Here we remove the unnecesary pieces by removing the extra + % coefficients. + coefs = zeros(size([coefX;coefY])-[8,0]); + coefs(1:2:end,:) = coefX(3:end-2,:); + coefs(2:2:end,:) = coefY(3:end-2,:); + + % Compute the actual arc length + seglen = sqrt(sum(diff(points,[],2).^2,1)); + cumarc = [0,cumsum(seglen)]; + + else + % The 1st derivatives at the end points and at the specified points + % are not equal, while analysed from both sides. Firstly compute + % the arclength of the point distribution. + seglen = sqrt(sum(diff(points,[],2).^2,1)); + cumarc = [0,cumsum(seglen)]; + + % Fit coefficients are obtained from the MATLAB(R) original pchip + % function, according to the desired derivative contiguity. + dp = find(d>0); + dpbig = find(diff(dp)>1); + dpbig = [dpbig,length(dp)]; + idx = dp(1):(dp(dpbig(1))+1); + temp = pchip(cumarc(idx),points(1,idx)); coefX = temp.coefs; + temp = pchip(cumarc(idx),points(2,idx)); coefY = temp.coefs; + for j=2:length(dpbig) + idx = dp(dpbig(j-1)+1):(dp(dpbig(j))+1); + temp = pchip(cumarc(idx),points(1,idx)); + coefX = [coefX;temp.coefs]; + temp = pchip(cumarc(idx),points(2,idx)); + coefY = [coefY;temp.coefs]; + end + + % Compiling the coefficients in a simple array. + coefs = zeros(size([coefX;coefY])); + coefs(1:2:end,:) = coefX(1:end,:); + coefs(2:2:end,:) = coefY(1:end,:); + + % Update the cumulative arclength + cumarc(:,d==0) = []; + end + + % Finally compute the piecewise polynomial. + curve = mkpp(cumarc,coefs,2); + + % Provide the differentiation array for later use. + diffarray = [3 0 0;0 2 0;0 0 1;0 0 0]; +end + +% If tq is given (or a print is required), compute the interpolation using +% the piecewise evaluation function provided in MATLAB(R) and then convert +% the parametrization into an arc-lenght one. +if tqgiven || print + step = (max(curve.breaks)-min(curve.breaks))/numel(tq)/30; + auxtq = min(curve.breaks):step:max(curve.breaks); + xyq = ppval(curve,auxtq); + tqp = pdearcl(auxtq,xyq,tq,0,1); + xyq = ppval(curve,tqp); +end + +% If the geometric parameters (perimeter and area) are required, compute +% them using some calculus. +if geomcalc + for ii = 1:curve.pieces + % Get the coefficients of the piecewise polynomial expresions of + % the parametric form. + coefX = curve.coefs(2*ii-1,:); + coefY = curve.coefs(2*ii,:); + + % Obtain the derivatives of the polynomials. + difX = coefX*diffarray; + difY = coefY*diffarray; + + % The length in the linear case is already computed, skip this bit. + if ~strcmpi(method,'linear') + % Define the function employed in the arc length and integrate + % it. + flen = @(t) sqrt(polyval(difX,t-curve.breaks(ii)).^2 ... + + polyval(difY,t-curve.breaks(ii)).^2); + len = len + integral(flen,curve.breaks(ii),curve.breaks(ii+1)); + end + % The area integral is computed here. + farea = @(t) polyval(conv(coefY,difX),... + t-curve.breaks(ii)); + area = area + integral(farea,curve.breaks(ii),... + curve.breaks(ii+1)); + + % The centroid is computed here + fcx = @(t) polyval(conv(coefX,conv(coefY,difX)),... + t-curve.breaks(ii)); + c(1) = c(1) + integral(fcx,curve.breaks(ii),... + curve.breaks(ii+1)); + fcy = @(t) polyval(conv(coefY,conv(coefX,difY)),... + t-curve.breaks(ii)); + c(2) = c(2) - integral(fcy,curve.breaks(ii),... + curve.breaks(ii+1)); + + % The area moments of inertia + fIxx = @(t) polyval(conv(coefY,conv(coefY,conv(coefY,difX))),... + t-curve.breaks(ii)); + Ixx = Ixx + integral(fIxx,curve.breaks(ii),... + curve.breaks(ii+1)); + fIyy = @(t) polyval(conv(coefX,conv(coefX,conv(coefX,difY))),... + t-curve.breaks(ii)); + Iyy = Iyy - integral(fIyy,curve.breaks(ii),... + curve.breaks(ii+1)); + fIxy = @(t) polyval(conv(coefY,conv(coefY,conv(coefX,difX))),... + t-curve.breaks(ii)); + Ixy = Ixy + integral(fIxy,curve.breaks(ii),... + curve.breaks(ii+1)); + end + c = c / area; + I = [1/3*Ixx,1/3*Iyy,1/2*Ixy]*sign(area); + area = abs(area); + + +end + +%% If required, print some figures to show what the algorithm did. +if print + figure + subplot(1,2,1) + plot(xyq(1,:),xyq(2,:),'*') + hold on + plot(points(1,:),points(2,:),'o') + plot(c(1),c(2),'x') + xlabel('x'), ylabel('y'), hold off, axis equal + title('Cartesian representation'), legend('Interpolation','Points',... + 'Centroid') + subplot(1,2,2) + plot(tq,xyq), hold on + for ii = (curve.breaks)/max(curve.breaks) + line([ii ii],ylim,'LineStyle','--','Color','k') + line([ii ii],ylim,'LineStyle','--','Color','k') + end + hold off, xlim([min(tq) max(tq)]), title('Parametric representation') + xlabel('t'), ylabel('x(t), y(t)'), legend('x(t)','y(t)') +end + +%% Process adequately the variables to be returned. +if nargout == 2 && ~tqgiven + varargout{1} = len; varargout{2} = area; +elseif nargout == 3 && ~tqgiven + varargout{1} = len; varargout{2} = area; varargout{3} = c; +elseif nargout == 3 + varargout{1} = xyq; + varargout{2} = len; varargout{3} = area; +elseif nargout == 4 && tqgiven + varargout{1} = xyq; + varargout{2} = len; varargout{3} = area; varargout{4} = c; +elseif nargout == 4 && ~tqgiven + varargout{1} = len; varargout{2} = area; varargout{3} = c; + varargout{4} = I; +elseif nargout == 5 + varargout{1} = xyq; + varargout{2} = len; varargout{3} = area; varargout{4} = c; + varargout{5} = I; +elseif nargout == 1 && pp + varargout{1} = curve; +else + varargout{1} = xyq; +end \ No newline at end of file diff --git a/lib/lofting.m b/lib/lofting.m new file mode 100644 index 0000000..93a2cb9 --- /dev/null +++ b/lib/lofting.m @@ -0,0 +1,30 @@ +function model = lofting(model,data_dir) +%Constructed Irregular 3D volumes from 2D contour curves +% example: +% Comsol_with_Matlab_Start; +% import com.comsol.model.util.* +% model = ModelUtil.create('Model1');% ModelUtil.remove('Model'); +% model.modelNode.create('mod1'); +% model.geom.create('geom1', 3); +% model.mesh.create('mesh1', 'geom1'); +% data_dir = pwd ; +% model = lofting(model,data_dir) + +CurveFiles = get_curveTxtFile(data_dir,[1:5],'curve'); +objectIPC = AndInterpolationCurve(model,CurveFiles); + +model.geom("geom1").create("loft1", "Loft"); +model.geom("geom1").feature("loft1").selection("profile").set(objectIPC); +model.geom("geom1").feature("loft1").set("facepartitioning", "grid"); + +% model.geom("geom1").create("pare1", "PartitionEdges"); + +try +model.component("mod1").geom("geom1").run(); +catch + warning('The automatic lofting failed, so the Partition Edges had to be added manually.'); +end + mphlaunch(model); + + +end \ No newline at end of file diff --git a/lib/plotSlice.m b/lib/plotSlice.m new file mode 100644 index 0000000..0e630ef --- /dev/null +++ b/lib/plotSlice.m @@ -0,0 +1,27 @@ + + +figure; +xyzID = ['X';'Y';'Z']; + +scatter(Pn(:,1),Pn( :,2),3,"filled"); +hold on; +plot(Pn(k,1),Pn(k,2),'g--','LineWidth',2); +hold on; +plot(PI(1,:),PI(2,:),'k','LineWidth',2); + + +title(['Silce' num2str(i)],'FontSize',12,'FontWeight','bold'); +xlabel([xyzID(planeID(1)) '(m)']); +ylabel([xyzID(planeID(2)) '(m)']); +l=legend('Point cloud slice','Point cloud boundary','Smooth boundary'); +set(l,'Box','off','FontSize',10); +set(gca,'color','none','linewidth',1,'FontSize',12,'FontWeight','bold'); +set(gcf,'Position', [713.8000 224.2000 404.8000 361.6000]); +box on; + +axis tight +axis equal; +% xlim([80,200]); +% ylim([-120,120]); +% set(gca,'color','none'); +% set(gcf,'color','none'); diff --git a/lib/smooth1q.m b/lib/smooth1q.m new file mode 100644 index 0000000..660fa34 --- /dev/null +++ b/lib/smooth1q.m @@ -0,0 +1,241 @@ +function [z,s] = smooth1q(y,s,varargin) + +%SMOOTH1Q Quick & easy smoothing. +% Z = SMOOTH1Q(Y,S) smoothes the data Y using a DCT- or FFT-based spline +% smoothing method. Non finite data (NaN or Inf) are treated as missing +% values. +% +% S is the smoothing parameter. It must be a real positive scalar. The +% larger S is, the smoother the output will be. If S is empty (i.e. S = +% []), it is automatically determined by minimizing the generalized +% cross-validation (GCV) score. +% +% Z = SMOOTH1Q(...,'robust') carries out a robust smoothing that +% minimizes the influence of outlying data. +% +% Z = SMOOTH1Q(...,'periodic') assumes that the data to be smoothed must +% be periodic. +% +% [Z,S] = SMOOTH1Q(...) also returns the calculated value for the +% smoothness parameter S so that you can fine-tune the smoothing +% subsequently if required. +% +% SMOOTH1Q is a simplified and quick version of SMOOTHN for 1-D data. If +% you want to smooth N-D arrays use SMOOTHN. +% +% Notes +% ----- +% 1) SMOOTH1Q works with regularly spaced data only. Use SMOOTH1 for non +% regularly spaced data. +% 2) The smoothness parameter used in this algorithm is determined +% automatically by minimizing the generalized cross-validation score. +% See the references for more details. +% +% References +% ---------- +% 1) Garcia D, Robust smoothing of gridded data in one and higher +% dimensions with missing values. Computational Statistics & Data +% Analysis, 2010. +% PDF download +% 2) Buckley MJ, Fast computation of a discretized thin-plate smoothing +% spline for image data. Biometrika, 1994. +% Link +% +% Examples: +% -------- +% % Simple curve +% x = linspace(0,100,200); +% y = cos(x/10)+(x/50).^2 + randn(size(x))/10; +% z = smooth1q(y,[]); +% plot(x,y,'r.',x,z,'k','LineWidth',2) +% axis tight square +% +% % Periodic curve with ouliers and missing data +% x = linspace(0,2*pi,300); +% y = cos(x)+ sin(2*x+1).^2 + randn(size(x))/5; +% y(150:155) = rand(1,6)*5; +% y(10:40) = NaN; +% subplot(121) +% z = smooth1q(y,1e3,'periodic'); +% plot(x,y,'r.',x,z,'k','LineWidth',2) +% axis tight square +% title('Non robust') +% subplot(122) +% z = smooth1q(y,1e3,'periodic','robust'); +% plot(x,y,'r.',x,z,'k','LineWidth',2) +% axis tight square +% title('Robust') +% +% % Limaon +% t = linspace(0,2*pi,300); +% x = cos(t).*(.5+cos(t)) + randn(size(t))*0.05; +% y = sin(t).*(.5+cos(t)) + randn(size(t))*0.05; +% z = smooth1q(complex(x,y),[],'periodic'); +% plot(x,y,'r.',real(z),imag(z),'k','linewidth',2) +% axis equal tight +% +% See also SMOOTHN, SMOOTH1. +% +% -- Damien Garcia -- 2012/08, revised 2014/02/26 +% website: www.BiomeCardio.com + +%-- Check input arguments +error(nargchk(2,4,nargin)); +assert(isvector(squeeze(y)),... + ['Y must be a 1-D array. Use SMOOTHN for non vector arrays.']) +if isempty(s) + isauto = 1; +else + assert(isnumeric(s),'S must be a numeric scalar') + assert(isscalar(s) && s>0,... + 'The smoothing parameter S must be a scalar >0') + isauto = 0; +end + +%-- Order (use m>=2, m = 2 is recommended) +m = 2; % Note: order of the smoothing process, can be modified + +%-- Options ('robust' and/or 'periodic') +isrobust = 0; method = 'dct'; % default options +%-- +if nargin>2 + assert(all(cellfun(@ischar,varargin)),... + 'The options must be ''robust'' and/or ''periodic''.') + varargin = lower(varargin); + if nargin==3 + idx = ismember({'robust','periodic'},varargin); + assert(any(idx),... + 'The options must be ''robust'' and/or ''periodic''.') + if idx(1), isrobust = 1; else method = 'fft'; end + else % nargin = 4 + assert(all(ismember(varargin,{'robust','periodic'})),... + 'The options must be ''robust'' and/or ''periodic''.') + isrobust = 1; + method = 'fft'; + end +end + +n = length(y); +siz0 = size(y); +y = y(:).'; + +%-- Weights +W0 = ones(siz0); +I = isfinite(y); % missing data (NaN or Inf values) +if any(~I) % replace the missing data (for faster convergence) + X = 1:n; + x = X(I); xi = X(~I); + y(~I) = interp1(x,y(I),xi,'linear','extrap'); +end +W0(~I) = 0; % weights for missing data are 0 +W = W0; + +%-- Eigenvalues +switch method + case 'dct' + Lambda = 2-2*cos((0:n-1)*pi/n); + case 'fft' + Lambda = 2-2*cos(2*(0:n-1)*pi/n); +end + +%-- Smoothing process +nr = 3; % Number of robustness iterations +for k = 0:nr*isrobust + if isrobust && k>0 + tmp = sqrt(1+16*s); + h = sqrt(1+tmp)/sqrt(2)/tmp; + W = W0.*bisquare(y,z,I,h); + end + if ~all(W==1) % then use an iterative method + tol = Inf; + zz = y; + while tol>1e-3 + switch method + case 'dct' + Y = dct(W.*(y-zz)+zz); + case 'fft' + Y = fft(W.*(y-zz)+zz); + end + if isauto + fminbnd(@GCVscore,-10,30,optimset('TolX',.1)); + else + Gamma = 1./(1+s*Lambda.^m); + switch method + case 'dct' + z = idct(Gamma.*Y); + case 'fft' + if isreal(y) + z = ifft(Gamma.*Y,'symmetric'); + else + z = ifft(Gamma.*Y); + end + end + end + tol = norm(zz-z)/norm(z); + zz = z; + end + + else %--- + % No missing values, non robust method => Direct fast method + %--- + switch method + case 'dct' + Y = dct(y); + case 'fft' + Y = fft(y); + end + if isauto + fminbnd(@GCVscore,-10,30,optimset('TolX',.1)); + else + Gamma = 1./(1+s*Lambda.^m); + end + switch method + case 'dct' + z = idct(Gamma.*Y); + case 'fft' + if isreal(y) + z = ifft(Gamma.*Y,'symmetric'); + else + z = ifft(Gamma.*Y); + end + end + end +end + +z = reshape(z,siz0); + + function GCVs = GCVscore(p) + s = 10^p; + Gamma = 1./(1+s*Lambda.^m); + if any(W) + switch method + case 'dct' + z = idct(Gamma.*Y); + case 'fft' + if isreal(y) + z = ifft(Gamma.*Y,'symmetric'); + else + z = ifft(Gamma.*Y); + end + end + RSS = norm(sqrt(W).*(y-z))^2; + else % No missing values, non robust method => Direct fast method + RSS = norm(Y.*(Gamma-1))^2; + end + TrH = sum(Gamma); + GCVs = RSS/(1-TrH/n)^2; + end + +end + +function W = bisquare(y,z,I,h) +r = y-z; % residuals +MAD = median(abs(r(I)-median(r(I)))); % median absolute deviation +u = abs(r/(1.4826*MAD)/sqrt(1-h)); % studentized residuals +W = (1-(u/4.685).^2).^2.*((u/4.685)<1); % bisquare weights +end \ No newline at end of file diff --git a/lib/smooth1qExample.m b/lib/smooth1qExample.m new file mode 100644 index 0000000..302e30d --- /dev/null +++ b/lib/smooth1qExample.m @@ -0,0 +1,30 @@ + % Simple curve + x = linspace(0,100,200); + y = cos(x/10)+(x/50).^2 + randn(size(x))/10; + z = smooth1q(y,[]); + plot(x,y,'r.',x,z,'k','LineWidth',2) + axis tight square + + % Periodic curve with ouliers and missing data + x = linspace(0,2*pi,300); + y = cos(x)+ sin(2*x+1).^2 + randn(size(x))/5; + y(150:155) = rand(1,6)*5; + y(10:40) = NaN; + subplot(121) + z = smooth1q(y,1e3,'periodic'); + plot(x,y,'r.',x,z,'k','LineWidth',2) + axis tight square + title('Non robust') + subplot(122) + z = smooth1q(y,1e3,'periodic','robust'); + plot(x,y,'r.',x,z,'k','LineWidth',2) + axis tight square + title('Robust') + + % Lima鏾n + t = linspace(0,2*pi,300); + x = cos(t).*(.5+cos(t)) + randn(size(t))*0.05; + y = sin(t).*(.5+cos(t)) + randn(size(t))*0.05; + z = smooth1q(complex(x,y),[],'periodic'); + plot(x,y,'r.',real(z),imag(z),'k','linewidth',2) +% axis equal tight \ No newline at end of file