Parallel Toolbox for data integration

Hi

The following fonction is used in a much larger code to extract strains from diffraction rings (from tiff images). It is the weekest point of the code, taking 98% of the computing time. I have also noticed that it only uses a single core, however we have 32-core machines. Jobs are running for days when they could be completed in minutes.

I have tried multple techniques, such as changing the loop in parfor, but I cannot seem to manage to spread the job onto multiple cores.

Is there anything I am not seing?

PS: the line starting with the "try" is the one taking all the computation time and that needs to be parallelized.

Many thanks

for azimuth = iSet.azimuth_min:iSet.azimuth_step:iSet.azimuth_max-iSet.azimuth_step
   
intensity_interpolated = zeros(round(((iSet.r_max-iSet.r_min)/iSet.r_step)+1),1);
    
    iSet.azimuth_sector_step =  0.1;
    iSet.azimuth_sector_min = azimuth+iSet.azimuth_sector_step/2;
    iSet.azimuth_sector_max = azimuth+iSet.azimuth_step-iSet.azimuth_sector_step/2;
    
    for azimuth_sector = iSet.azimuth_sector_min:iSet.azimuth_sector_step:iSet.azimuth_sector_max% degrees
        
        intensity_temp_interpolated = zeros(round(((iSet.r_max-iSet.r_min)/iSet.r_step)+1),1);
        
        for r = iSet.r_min:iSet.r_step:iSet.r_max % in pixels
            
            r_index = round(((r-iSet.r_min)/iSet.r_step)+1);
            
            x(r_index) = r.*cos(azimuth_sector*pi./180) +  XCentreGuess;
            y(r_index) = r.*sin(azimuth_sector*pi./180) +  YCentreGuess;
            
            floor_x = floor(x(r_index));  floor_y = floor(y(r_index));
            ceil_x = ceil(x(r_index)); ceil_y = ceil(y(r_index));
            
            % using biliear interpolation (http://en.wikipedia.org/wiki/Bilinear_interpolation)
            try
                intensity_temp_interpolated(r_index) = area_pattern(floor_y,floor_x).*(1-(x(r_index)-floor_x)).*(1-(y(r_index)-floor_y)) ...
                    + area_pattern(ceil_y,floor_x).*(x(r_index)-floor_x).*(1-(y(r_index)-floor_y)) ...
                    + area_pattern(floor_y,ceil_x).*(1-(x(r_index)-floor_x)).*(y(r_index)-floor_y) ...
                    + area_pattern(ceil_y,ceil_x).*(x(r_index)-floor_x).*(y(r_index)-floor_y);
            catch
                intensity_temp_interpolated(r_index) = NaN;
            end
                        
        end
        
        intensity_interpolated = intensity_temp_interpolated + intensity_interpolated;
       
    end
            % draw the sector that was calculated and superimpose to the 2D diffraction
            
            if iSet.DrawFittedSector == 1
                        fDrawFittedSector(area_pattern,iSet.r_min,iSet.r_step,iSet.r_max,iSet.azimuth_sector_min,iSet.azimuth_sector_max,XCentreGuess,YCentreGuess)
            end
            
            azimuth_index = round(((azimuth-iSet.azimuth_min)/iSet.azimuth_step)+1);
            
    intensity_sector_average(:,azimuth_index) =  intensity_interpolated./(((iSet.azimuth_sector_max-iSet.azimuth_sector_min)/iSet.azimuth_sector_step) + 1);
end %%% new %%%


      
end