1 subroutine parsur(iopt,ipar,idim,mu,u,mv,v,f,s,nuest,nvest, 2 * nu,tu,nv,tv,c,fp,wrk,lwrk,iwrk,kwrk,ier) 3c given the set of ordered points f(i,j) in the idim-dimensional space, 4c corresponding to grid values (u(i),v(j)) ,i=1,...,mu ; j=1,...,mv, 5c parsur determines a smooth approximating spline surface s(u,v) , i.e. 6c f1 = s1(u,v) 7c ... u(1) <= u <= u(mu) ; v(1) <= v <= v(mv) 8c fidim = sidim(u,v) 9c with sl(u,v), l=1,2,...,idim bicubic spline functions with common 10c knots tu(i),i=1,...,nu in the u-variable and tv(j),j=1,...,nv in the 11c v-variable. 12c in addition, these splines will be periodic in the variable u if 13c ipar(1) = 1 and periodic in the variable v if ipar(2) = 1. 14c if iopt=-1, parsur determines the least-squares bicubic spline 15c surface according to a given set of knots. 16c if iopt>=0, the number of knots of s(u,v) and their position 17c is chosen automatically by the routine. the smoothness of s(u,v) is 18c achieved by minimalizing the discontinuity jumps of the derivatives 19c of the splines at the knots. the amount of smoothness of s(u,v) is 20c determined by the condition that 21c fp=sumi=1,mu(sumj=1,mv(dist(f(i,j)-s(u(i),v(j)))**2))<=s, 22c with s a given non-negative constant. 23c the fit s(u,v) is given in its b-spline representation and can be 24c evaluated by means of routine surev. 25c 26c calling sequence: 27c call parsur(iopt,ipar,idim,mu,u,mv,v,f,s,nuest,nvest,nu,tu, 28c * nv,tv,c,fp,wrk,lwrk,iwrk,kwrk,ier) 29c 30c parameters: 31c iopt : integer flag. unchanged on exit. 32c on entry iopt must specify whether a least-squares surface 33c (iopt=-1) or a smoothing surface (iopt=0 or 1)must be 34c determined. 35c if iopt=0 the routine will start with the initial set of 36c knots needed for determining the least-squares polynomial 37c surface. 38c if iopt=1 the routine will continue with the set of knots 39c found at the last call of the routine. 40c attention: a call with iopt=1 must always be immediately 41c preceded by another call with iopt = 1 or iopt = 0. 42c ipar : integer array of dimension 2. unchanged on exit. 43c on entry ipar(1) must specify whether (ipar(1)=1) or not 44c (ipar(1)=0) the splines must be periodic in the variable u. 45c on entry ipar(2) must specify whether (ipar(2)=1) or not 46c (ipar(2)=0) the splines must be periodic in the variable v. 47c idim : integer. on entry idim must specify the dimension of the 48c surface. 1 <= idim <= 3. unchanged on exit. 49c mu : integer. on entry mu must specify the number of grid points 50c along the u-axis. unchanged on exit. 51c mu >= mumin where mumin=4-2*ipar(1) 52c u : real array of dimension at least (mu). before entry, u(i) 53c must be set to the u-co-ordinate of the i-th grid point 54c along the u-axis, for i=1,2,...,mu. these values must be 55c supplied in strictly ascending order. unchanged on exit. 56c mv : integer. on entry mv must specify the number of grid points 57c along the v-axis. unchanged on exit. 58c mv >= mvmin where mvmin=4-2*ipar(2) 59c v : real array of dimension at least (mv). before entry, v(j) 60c must be set to the v-co-ordinate of the j-th grid point 61c along the v-axis, for j=1,2,...,mv. these values must be 62c supplied in strictly ascending order. unchanged on exit. 63c f : real array of dimension at least (mu*mv*idim). 64c before entry, f(mu*mv*(l-1)+mv*(i-1)+j) must be set to the 65c l-th co-ordinate of the data point corresponding to the 66c the grid point (u(i),v(j)) for l=1,...,idim ,i=1,...,mu 67c and j=1,...,mv. unchanged on exit. 68c if ipar(1)=1 it is expected that f(mu*mv*(l-1)+mv*(mu-1)+j) 69c = f(mu*mv*(l-1)+j), l=1,...,idim ; j=1,...,mv 70c if ipar(2)=1 it is expected that f(mu*mv*(l-1)+mv*(i-1)+mv) 71c = f(mu*mv*(l-1)+mv*(i-1)+1), l=1,...,idim ; i=1,...,mu 72c s : real. on entry (if iopt>=0) s must specify the smoothing 73c factor. s >=0. unchanged on exit. 74c for advice on the choice of s see further comments 75c nuest : integer. unchanged on exit. 76c nvest : integer. unchanged on exit. 77c on entry, nuest and nvest must specify an upper bound for the 78c number of knots required in the u- and v-directions respect. 79c these numbers will also determine the storage space needed by 80c the routine. nuest >= 8, nvest >= 8. 81c in most practical situation nuest = mu/2, nvest=mv/2, will 82c be sufficient. always large enough are nuest=mu+4+2*ipar(1), 83c nvest = mv+4+2*ipar(2), the number of knots needed for 84c interpolation (s=0). see also further comments. 85c nu : integer. 86c unless ier=10 (in case iopt>=0), nu will contain the total 87c number of knots with respect to the u-variable, of the spline 88c surface returned. if the computation mode iopt=1 is used, 89c the value of nu should be left unchanged between subsequent 90c calls. in case iopt=-1, the value of nu should be specified 91c on entry. 92c tu : real array of dimension at least (nuest). 93c on successful exit, this array will contain the knots of the 94c splines with respect to the u-variable, i.e. the position of 95c the interior knots tu(5),...,tu(nu-4) as well as the position 96c of the additional knots tu(1),...,tu(4) and tu(nu-3),..., 97c tu(nu) needed for the b-spline representation. 98c if the computation mode iopt=1 is used,the values of tu(1) 99c ...,tu(nu) should be left unchanged between subsequent calls. 100c if the computation mode iopt=-1 is used, the values tu(5), 101c ...tu(nu-4) must be supplied by the user, before entry. 102c see also the restrictions (ier=10). 103c nv : integer. 104c unless ier=10 (in case iopt>=0), nv will contain the total 105c number of knots with respect to the v-variable, of the spline 106c surface returned. if the computation mode iopt=1 is used, 107c the value of nv should be left unchanged between subsequent 108c calls. in case iopt=-1, the value of nv should be specified 109c on entry. 110c tv : real array of dimension at least (nvest). 111c on successful exit, this array will contain the knots of the 112c splines with respect to the v-variable, i.e. the position of 113c the interior knots tv(5),...,tv(nv-4) as well as the position 114c of the additional knots tv(1),...,tv(4) and tv(nv-3),..., 115c tv(nv) needed for the b-spline representation. 116c if the computation mode iopt=1 is used,the values of tv(1) 117c ...,tv(nv) should be left unchanged between subsequent calls. 118c if the computation mode iopt=-1 is used, the values tv(5), 119c ...tv(nv-4) must be supplied by the user, before entry. 120c see also the restrictions (ier=10). 121c c : real array of dimension at least (nuest-4)*(nvest-4)*idim. 122c on successful exit, c contains the coefficients of the spline 123c approximation s(u,v) 124c fp : real. unless ier=10, fp contains the sum of squared 125c residuals of the spline surface returned. 126c wrk : real array of dimension (lwrk). used as workspace. 127c if the computation mode iopt=1 is used the values of 128c wrk(1),...,wrk(4) should be left unchanged between subsequent 129c calls. 130c lwrk : integer. on entry lwrk must specify the actual dimension of 131c the array wrk as declared in the calling (sub)program. 132c lwrk must not be too small. 133c lwrk >= 4+nuest*(mv*idim+11+4*ipar(1))+nvest*(11+4*ipar(2))+ 134c 4*(mu+mv)+q*idim where q is the larger of mv and nuest. 135c iwrk : integer array of dimension (kwrk). used as workspace. 136c if the computation mode iopt=1 is used the values of 137c iwrk(1),.,iwrk(3) should be left unchanged between subsequent 138c calls. 139c kwrk : integer. on entry kwrk must specify the actual dimension of 140c the array iwrk as declared in the calling (sub)program. 141c kwrk >= 3+mu+mv+nuest+nvest. 142c ier : integer. unless the routine detects an error, ier contains a 143c non-positive value on exit, i.e. 144c ier=0 : normal return. the surface returned has a residual sum of 145c squares fp such that abs(fp-s)/s <= tol with tol a relat- 146c ive tolerance set to 0.001 by the program. 147c ier=-1 : normal return. the spline surface returned is an 148c interpolating surface (fp=0). 149c ier=-2 : normal return. the surface returned is the least-squares 150c polynomial surface. in this extreme case fp gives the 151c upper bound for the smoothing factor s. 152c ier=1 : error. the required storage space exceeds the available 153c storage space, as specified by the parameters nuest and 154c nvest. 155c probably causes : nuest or nvest too small. if these param- 156c eters are already large, it may also indicate that s is 157c too small 158c the approximation returned is the least-squares surface 159c according to the current set of knots. the parameter fp 160c gives the corresponding sum of squared residuals (fp>s). 161c ier=2 : error. a theoretically impossible result was found during 162c the iteration process for finding a smoothing surface with 163c fp = s. probably causes : s too small. 164c there is an approximation returned but the corresponding 165c sum of squared residuals does not satisfy the condition 166c abs(fp-s)/s < tol. 167c ier=3 : error. the maximal number of iterations maxit (set to 20 168c by the program) allowed for finding a smoothing surface 169c with fp=s has been reached. probably causes : s too small 170c there is an approximation returned but the corresponding 171c sum of squared residuals does not satisfy the condition 172c abs(fp-s)/s < tol. 173c ier=10 : error. on entry, the input data are controlled on validity 174c the following restrictions must be satisfied. 175c -1<=iopt<=1, 0<=ipar(1)<=1, 0<=ipar(2)<=1, 1 <=idim<=3 176c mu >= 4-2*ipar(1),mv >= 4-2*ipar(2), nuest >=8, nvest >= 8, 177c kwrk>=3+mu+mv+nuest+nvest, 178c lwrk >= 4+nuest*(mv*idim+11+4*ipar(1))+nvest*(11+4*ipar(2)) 179c +4*(mu+mv)+max(nuest,mv)*idim 180c u(i-1)<u(i),i=2,..,mu, v(i-1)<v(i),i=2,...,mv 181c if iopt=-1: 8<=nu<=min(nuest,mu+4+2*ipar(1)) 182c u(1)<tu(5)<tu(6)<...<tu(nu-4)<u(mu) 183c 8<=nv<=min(nvest,mv+4+2*ipar(2)) 184c v(1)<tv(5)<tv(6)<...<tv(nv-4)<v(mv) 185c the schoenberg-whitney conditions, i.e. there must 186c be subset of grid co-ordinates uu(p) and vv(q) such 187c that tu(p) < uu(p) < tu(p+4) ,p=1,...,nu-4 188c tv(q) < vv(q) < tv(q+4) ,q=1,...,nv-4 189c (see fpchec or fpchep) 190c if iopt>=0: s>=0 191c if s=0: nuest>=mu+4+2*ipar(1) 192c nvest>=mv+4+2*ipar(2) 193c if one of these conditions is found to be violated,control 194c is immediately repassed to the calling program. in that 195c case there is no approximation returned. 196c 197c further comments: 198c by means of the parameter s, the user can control the tradeoff 199c between closeness of fit and smoothness of fit of the approximation. 200c if s is too large, the surface will be too smooth and signal will be 201c lost ; if s is too small the surface will pick up too much noise. in 202c the extreme cases the program will return an interpolating surface 203c if s=0 and the constrained least-squares polynomial surface if s is 204c very large. between these extremes, a properly chosen s will result 205c in a good compromise between closeness of fit and smoothness of fit. 206c to decide whether an approximation, corresponding to a certain s is 207c satisfactory the user is highly recommended to inspect the fits 208c graphically. 209c recommended values for s depend on the accuracy of the data values. 210c if the user has an idea of the statistical errors on the data, he 211c can also find a proper estimate for s. for, by assuming that, if he 212c specifies the right s, parsur will return a surface s(u,v) which 213c exactly reproduces the surface underlying the data he can evaluate 214c the sum(dist(f(i,j)-s(u(i),v(j)))**2) to find a good estimate for s. 215c for example, if he knows that the statistical errors on his f(i,j)- 216c values is not greater than 0.1, he may expect that a good s should 217c have a value not larger than mu*mv*(0.1)**2. 218c if nothing is known about the statistical error in f(i,j), s must 219c be determined by trial and error, taking account of the comments 220c above. the best is then to start with a very large value of s (to 221c determine the le-sq polynomial surface and the corresponding upper 222c bound fp0 for s) and then to progressively decrease the value of s 223c ( say by a factor 10 in the beginning, i.e. s=fp0/10,fp0/100,... 224c and more carefully as the approximation shows more detail) to 225c obtain closer fits. 226c to economize the search for a good s-value the program provides with 227c different modes of computation. at the first call of the routine, or 228c whenever he wants to restart with the initial set of knots the user 229c must set iopt=0. 230c if iopt = 1 the program will continue with the knots found at 231c the last call of the routine. this will save a lot of computation 232c time if parsur is called repeatedly for different values of s. 233c the number of knots of the surface returned and their location will 234c depend on the value of s and on the complexity of the shape of the 235c surface underlying the data. if the computation mode iopt = 1 236c is used, the knots returned may also depend on the s-values at 237c previous calls (if these were smaller). therefore, if after a number 238c of trials with different s-values and iopt=1,the user can finally 239c accept a fit as satisfactory, it may be worthwhile for him to call 240c parsur once more with the chosen value for s but now with iopt=0. 241c indeed, parsur may then return an approximation of the same quality 242c of fit but with fewer knots and therefore better if data reduction 243c is also an important objective for the user. 244c the number of knots may also depend on the upper bounds nuest and 245c nvest. indeed, if at a certain stage in parsur the number of knots 246c in one direction (say nu) has reached the value of its upper bound 247c (nuest), then from that moment on all subsequent knots are added 248c in the other (v) direction. this may indicate that the value of 249c nuest is too small. on the other hand, it gives the user the option 250c of limiting the number of knots the routine locates in any direction 251c for example, by setting nuest=8 (the lowest allowable value for 252c nuest), the user can indicate that he wants an approximation with 253c splines which are simple cubic polynomials in the variable u. 254c 255c other subroutines required: 256c fppasu,fpchec,fpchep,fpknot,fprati,fpgrpa,fptrnp,fpback, 257c fpbacp,fpbspl,fptrpe,fpdisc,fpgivs,fprota 258c 259c author: 260c p.dierckx 261c dept. computer science, k.u. leuven 262c celestijnenlaan 200a, b-3001 heverlee, belgium. 263c e-mail : Paul.Dierckx@cs.kuleuven.ac.be 264c 265c latest update : march 1989 266c 267c .. 268c ..scalar arguments.. 269 real*8 s,fp 270 integer iopt,idim,mu,mv,nuest,nvest,nu,nv,lwrk,kwrk,ier 271c ..array arguments.. 272 real*8 u(mu),v(mv),f(mu*mv*idim),tu(nuest),tv(nvest), 273 * c((nuest-4)*(nvest-4)*idim),wrk(lwrk) 274 integer ipar(2),iwrk(kwrk) 275c ..local scalars.. 276 real*8 tol,ub,ue,vb,ve,peru,perv 277 integer i,j,jwrk,kndu,kndv,knru,knrv,kwest,l1,l2,l3,l4, 278 * lfpu,lfpv,lwest,lww,maxit,nc,mf,mumin,mvmin 279c ..function references.. 280 integer max0 281c ..subroutine references.. 282c fppasu,fpchec,fpchep 283c .. 284c we set up the parameters tol and maxit. 285 maxit = 20 286 tol = 0.1e-02 287c before starting computations a data check is made. if the input data 288c are invalid, control is immediately repassed to the calling program. 289 ier = 10 290 if(iopt.lt.(-1) .or. iopt.gt.1) go to 200 291 if(ipar(1).lt.0 .or. ipar(1).gt.1) go to 200 292 if(ipar(2).lt.0 .or. ipar(2).gt.1) go to 200 293 if(idim.le.0 .or. idim.gt.3) go to 200 294 mumin = 4-2*ipar(1) 295 if(mu.lt.mumin .or. nuest.lt.8) go to 200 296 mvmin = 4-2*ipar(2) 297 if(mv.lt.mvmin .or. nvest.lt.8) go to 200 298 mf = mu*mv 299 nc = (nuest-4)*(nvest-4) 300 lwest = 4+nuest*(mv*idim+11+4*ipar(1))+nvest*(11+4*ipar(2))+ 301 * 4*(mu+mv)+max0(nuest,mv)*idim 302 kwest = 3+mu+mv+nuest+nvest 303 if(lwrk.lt.lwest .or. kwrk.lt.kwest) go to 200 304 do 10 i=2,mu 305 if(u(i-1).ge.u(i)) go to 200 306 10 continue 307 do 20 i=2,mv 308 if(v(i-1).ge.v(i)) go to 200 309 20 continue 310 if(iopt.ge.0) go to 100 311 if(nu.lt.8 .or. nu.gt.nuest) go to 200 312 ub = u(1) 313 ue = u(mu) 314 if (ipar(1).ne.0) go to 40 315 j = nu 316 do 30 i=1,4 317 tu(i) = ub 318 tu(j) = ue 319 j = j-1 320 30 continue 321 call fpchec(u,mu,tu,nu,3,ier) 322 if(ier.ne.0) go to 200 323 go to 60 324 40 l1 = 4 325 l2 = l1 326 l3 = nu-3 327 l4 = l3 328 peru = ue-ub 329 tu(l2) = ub 330 tu(l3) = ue 331 do 50 j=1,3 332 l1 = l1+1 333 l2 = l2-1 334 l3 = l3+1 335 l4 = l4-1 336 tu(l2) = tu(l4)-peru 337 tu(l3) = tu(l1)+peru 338 50 continue 339 call fpchep(u,mu,tu,nu,3,ier) 340 if(ier.ne.0) go to 200 341 60 if(nv.lt.8 .or. nv.gt.nvest) go to 200 342 vb = v(1) 343 ve = v(mv) 344 if (ipar(2).ne.0) go to 80 345 j = nv 346 do 70 i=1,4 347 tv(i) = vb 348 tv(j) = ve 349 j = j-1 350 70 continue 351 call fpchec(v,mv,tv,nv,3,ier) 352 if(ier.ne.0) go to 200 353 go to 150 354 80 l1 = 4 355 l2 = l1 356 l3 = nv-3 357 l4 = l3 358 perv = ve-vb 359 tv(l2) = vb 360 tv(l3) = ve 361 do 90 j=1,3 362 l1 = l1+1 363 l2 = l2-1 364 l3 = l3+1 365 l4 = l4-1 366 tv(l2) = tv(l4)-perv 367 tv(l3) = tv(l1)+perv 368 90 continue 369 call fpchep(v,mv,tv,nv,3,ier) 370 if (ier.eq.0) go to 150 371 go to 200 372 100 if(s.lt.0.) go to 200 373 if(s.eq.0. .and. (nuest.lt.(mu+4+2*ipar(1)) .or. 374 * nvest.lt.(mv+4+2*ipar(2))) )go to 200 375 ier = 0 376c we partition the working space and determine the spline approximation 377 150 lfpu = 5 378 lfpv = lfpu+nuest 379 lww = lfpv+nvest 380 jwrk = lwrk-4-nuest-nvest 381 knru = 4 382 knrv = knru+mu 383 kndu = knrv+mv 384 kndv = kndu+nuest 385 call fppasu(iopt,ipar,idim,u,mu,v,mv,f,mf,s,nuest,nvest, 386 * tol,maxit,nc,nu,tu,nv,tv,c,fp,wrk(1),wrk(2),wrk(3),wrk(4), 387 * wrk(lfpu),wrk(lfpv),iwrk(1),iwrk(2),iwrk(3),iwrk(knru), 388 * iwrk(knrv),iwrk(kndu),iwrk(kndv),wrk(lww),jwrk,ier) 389 200 return 390 end 391 392