Clear Filters
Clear Filters

N-Dimensional smoothing spline weighting

11 views (last 30 days)
Avery Berman
Avery Berman on 28 Mar 2016
Answered: Avery Berman on 11 Dec 2018
Is anybody aware of how to use the weighting option of the csaps function for N-dimensional matrices? As an example, I can currently call csaps to fit a smoothing spline to a 2D MxN image, I, with a command like:
x = {1:M, 1:N};
spl_smooth = csaps(x, I);
I would like to fit a smoothing spline using a weighting matrix that has the same dimensions as the input image but it does not appear that I can do this with csaps since the weighting input argument must be a cell array with the same number of entries as the first input, x. This makes the size of the weighting argument {1xM, 1xN}, whereas I want it to be MxN.
Thank you for any help you can offer,

Answers (2)

Unai San Miguel
Unai San Miguel on 10 Dec 2018
You can use a cell of two set of weights, one in each direction of the input
csaps({1:M, 1:N}, I, [], {wi, wj});
with wi an array of weights of size [1, M] and wj an array of weights of size [1, N]. In this way at least you have constant-slices weighting.
Avery Berman
Avery Berman on 11 Dec 2018
Thanks again for the suggestion Unai. I believe that is what I originally did as well. The issue with masking the matrix itself prior to applying the weighting, is that the spline will now attempt to fit the untrusted regions to 0s. This has the result of artifactually decreasing the smoothed matrix in the regions that were trusted. I was adequately satisfied with my workaround, although it would have been nice if I could have done that all in one step with an element-by-element weighting array.
Unai San Miguel
Unai San Miguel on 11 Dec 2018
And varying the p parameters from {0.9, 0.9} as defaults to some other values?
A smoothing spline will always decrease the smoothed matrix in the trusted regions (and increase it in the untrusted ones), won't it?
I would try to reason first in 1D. If you could make up an example in 1D, that is, 1-valued, 1-variate functions maybe I could understand better your case. Not meaning I would be able to help, though...

Sign in to comment.

Avery Berman
Avery Berman on 11 Dec 2018
Hi Unai,
Sorry for the confusion, my original posting was a couple of years ago, so I was not entirely clear on the details myself! You are totally right that the smoothing should decrease the signal if the untrusted regions naturally have low signal. What I was missing was that I was dealing with complex data with very slowly varying phase. In some regions, the signal magnitude would drop to very low levels, in which case the phase data is very noisy. My intention was to recapture the phase signal, not the magnitude. I was using the magnitude to define my mask of pixels (or array elements) to trust, and then interpolating the phase in the non-masked regions, then smoothing the phase.
I've included a simple 1D example of this where I have a true signal of length 50 whose magnitude is either 0 on the edges and 1 throughout the middle (from 15 to 40), except a small region where the intensity drops down to 0.05 (from 20 to 24). Throughout the middle, the phase is linearly increasing within the central region. Finally, I've added Gaussian distributed noise to the signal. My intention is to recover the phase in the region where the signal magnitude was low, which made the phase very noisy. If I fit the spline in the central region, then the smoothed spline dips down where it shouldn't**. If I first interpolate the elements with I < 0.10 using griddedInterpolant (the 1D equivalent of scatteredInterpolant), the interpolant helps fill in the central region, then I could smooth this with a spline if I wanted to. Hopefully, however, this illustrates how I cannot just set the phase to 0 in the low magnitude region since the spline will attempt to fit this area.
**Note that in the 1D case here, I could have used a weighting input argument to csaps but that is not possible in my 2D (or 3D) applications, so for the sake of the example, I've ignored the weighting.
Thanks again for your feedback,

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!