Optimization of FIT function in a FOR loop

My current project work is based on processing tiff images. I have x images each of size m x n pixels.

I am using Matlab FIT function to fit a customised equation using data from all x images at each pixel (mxn). From the fit, three parameters need to be calculated.

I am using 2 for loop to calculate the fit for each pixel. The program is extremely slow (> 6 hours) during execution of the 2 for loops even with parallel processing. Preallocation is already included in the code.

Is there a way to call the FIT function just once and it executes the fitting for all pixels?

The relevant code is attached here, fitt() defines all the options for fit V1 & S1 are input parameters, X1 Y1 is dimensions of the images.

parfor i = 1:X1
        [ft1, opts1] = fitt();
       for j = 1:Y1

              [result1, gof1] = fit(V1(:,i,j),S1(:,i,j),ft1,opts1);
              A1(i,j)=result1.C3;
              A2(i,j)=result1.C1;
              A3(i,j)=result1.C2;

          end
      end

I am using Matlab version R2012a on Windows Vista OS.

Thanks

Answer 1

[Sean de Wolski]

A few questions and thoughts:

• What are you actually fitting and how big are your images?
• Also, have you tried LSQCURVEFIT in the Optimization Toolbox? I think you might see a speedup with it.
• Another thing to consider is how much variation do you expect in the coefficients of adjacent or near adjacent pixels? If the variation is not expected to be much, use the results from the previous iteration as your initial guess in the next iteration. This will hopefully give the optimizer a good starting point so that it can converge with fewer iterations.
• Do you really need to do this at each pixel? Downsampling by a factor of two in each dimension will get you a 4x speedup pretty effortlessly.

[Pooja C]

The fit equation is sum of three variables - two exponentials and one constant. The original image is 1K x1K but I have cropped and downsampled the image to 350x350 pixels.

The coefficients are specific to material quality of the sample imaged. Poor or very good quality regions may have large variations.I can try using results from adjacent pixels as a starting point.

I can try LSQCURVEFIT function. Is this faster than FIT function? Any other advantages?

Thanks so much for your suggestions.

[Sean de Wolski]

I haven't timed it but expect it to be faster as it is just a solver, it does not compute other things as well.

[Pooja C]

Sorry for my late reply.

The LSQCURVEFIT function is definitely faster and reduces the computation time by almost 50%. I checked the values computed with FIT and LSQCURVEFIT and they are very close in terms of min, max and mean and also pixelwise subtraction or division . Thanks so much for that tip.

I had some problem using the start value of the three parameters from the values of the adjacent pixels. I am using parallel processing and it is definitely not feasible to use there in straightforward way.

[Sean de Wolski]

Hi Pooja,

Glad to hear lsqcurvefit sped it up. You are certainly right that PARFOR requires independent iterations. However, there are approaches you could take to get around this. Perhaps a parfor loop with a regular for loop nested inside of it. You could also do a coarse search, say look at every 3rd or 4th pixel in each dimension. Save these results, and then use that as a starting point for all of the neighbor pixels in a second round of parfor. This would alleviate the dependency between iterations.

[Matt J]

It might help further to use LSQCURVEFIT with a user-supplied Jacobian calculation, if this is not being done already,

options=optimoptions(@lsqcurvefit,'Jacobian','on')

It sounds like this could cut down on the number of exp() operations by at least a few factors, and might also accelerate convergence.

[Pooja C]

Hi Sean, I will try your suggestion soon. I am working on least number of images needed for curve fitting and LSQCURVEFIT has indeed made the computation so much faster.

Hi Matt J, The Jacobian showed little improvement in time reduction.

Thanks to both of you.

[Matt J]

Hi Matt J, The Jacobian showed little improvement in time reduction.

It's more than just turning Jacobian on. There is code optimization to be done. You have to recycle the exponential expressions used in the objective function calculation when computing the Jacobian.

[Pooja C]

Thanks Matt. Will try it.