N-Dimensional smoothing spline weighting

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,

Avery

Answer 1

[Unai San Miguel]

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]

Thanks Unai!

The problem with this is that it still does not provide independent control of each pixel's weighting. What I was mostly interested in was being able to smooth an image and to only calculate the smoothing function using a mask of pixels that I trusted. My workaround was to first use scatteredInterpolant to fill in the pixels that I did not trust (by linear interpolation), then I followed that with an unweighted spline fit to effectively smooth the entire image.

For example, if I only trusted pixels whose intensity was greater than 5% of the max intensity:

mask_in = ( I > (0.05 * max(I(:))) );  % I is the input image
x = 1:size(I,1);
y = 1:size(I,2);
[xx, yy] = meshgrid(x,y);

% Create a scattered interpolant using mask_in and apply to all xx, yy to interpolate
ScattInt = scatteredInterpolant( xx(mask_in), yy(mask_in), I(mask_in), 'linear' );
I_int = ScattInt(xx, yy);

% Fit a smoothing spline to the interpolated image
[spl_fit, p] = csaps( {x, y}, I_int );
% Apply the spline
I_smooth  = fnval(spl_fit, {x, y});

I cannot guarantee that this will run without error as this is an adaptation of my actual code, but hopefully this at least conveys the message.

Avery

[Unai San Miguel]

I'm sorry but I don't usually work with images, and I am usually very confused with meshgrid (I prefer ndgrid...)

If you want to weight more some regions of your image you could directly multiply those regions by a weighting function,as in

mask_in = ( I > (0.5 * max(I(:))) );  % I is the input image
x = 1:size(I,1);
y = 1:size(I,2);

[spl_fit2, p2] = csaps({x, y}, I .* mask_in)
I_smooth2 = fnval(spl_fit2, {x, y});

Is this something you could use?

[Avery Berman]

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]

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...

Answer 2

[Avery Berman]

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,

Avery