Finding the points inside an ellipse, a utility function to be used with matlab regionprops

% borrowed the idea from Steve's blog 
% and wrote the function I needed, here is a link to Steve's blog

url = 'http://blogs.mathworks.com/images/steve/2010/rice_binary.png';
bw = imread(url);

dataArray = zeros(size(bw));
dataArray = rand(size(bw));
s = regionprops(bw, 'Orientation', 'MajorAxisLength', ...
    'MinorAxisLength', 'Eccentricity', 'Centroid');
imagesc(bw)
hold on

phi = [0,pi/2, pi, 2*pi*3/4.0];
cosphi = cos(phi);
sinphi = sin(phi);
for k = 1:length(s)
    xbar = s(k).Centroid(1);
    ybar = s(k).Centroid(2);

    a = s(k).MajorAxisLength/2;
    b = s(k).MinorAxisLength/2;

    theta = pi*s(k).Orientation/180;   
    R = [ cos(theta)   sin(theta)
         -sin(theta)   cos(theta)];
    xy = [a*cosphi; b*sinphi];
    xy = R*xy;
    x = xy(1,:) + xbar;
    y = xy(2,:) + ybar;
    plot(x,y,'*g'); 
    [avg, totalElem, jx, jy] = avgInEllipse(xbar, ybar, a, b,...
        s(k).Orientation, bw, dataArray);
    xvec = zeros(jx.size(),1);
    yvec = zeros(jy.size(),1);
    for tt = 0:jx.size()-1
        xvec(tt+1) = jx.elementAt(tt);
        yvec(tt+1) = jy.elementAt(tt);
    end
    plot(xvec, yvec, '.m')
    disp(['total no of interior points :' num2str(totalElem)]);
    disp(['Average :' num2str(avg)]);     
    
end



% this part is directly copied from Steve's blog, I just inserted comments

% create 50 points around the circumference of the unit circle
phi = linspace(0,2*pi,50);
cosphi = cos(phi);
sinphi = sin(phi);

for k = 1:length(s)
    xbar = s(k).Centroid(1);
    ybar = s(k).Centroid(2);
    % get radii of the two circles centered at (0,0)
    a = s(k).MajorAxisLength/2;
    b = s(k).MinorAxisLength/2;
    % convert Orientation angle to radian value
    theta = pi*s(k).Orientation/180;
    
    % rotation matrix
    R = [ cos(theta)   sin(theta)
         -sin(theta)   cos(theta)];
    % a*cosphi gives points on the x-axis of a circle
    % drawn around (0,0) point with a radius of 'a'
    % b*sinphi gives points on the y-axis of a circle
    % drawn around (0,0) point with a radius of 'b'    
    xy = [a*cosphi; b*sinphi];
    % rotate the points to match orientation
    xy = R*xy;
    % translate the points that were generated around (0,0) point
    % to the new center (xbar, ybar) aka centroid
    x = xy(1,:) + xbar;
    y = xy(2,:) + ybar;

    plot(x,y,'r','LineWidth',2);    
end

hold off

Input Image

output image, enclosed points drawn in purple

output image, two elements are zoomed in

function [avg, totalElem, jvector_x, jvector_y] = avgInEllipse(cx, cy, a, b,...
    orientation, bw, dataArray)
    % Author : Mohammad A. Khan
    % =======================  inputs ==========================
    % this function is to be used in conjunction with regionprops
    % cx, cy are the x and y coordinated of the centroid
    % a, b are the 1/2 lengths of the major and minor axes of the 
    % ellipse returned by regionprops, orientation is the orientation
    % returned by regionprops, which is the degree of angle between
    % ellipse's major axis and x-axis, bw is the binary array where
    % elements were detected, it can also be though of the mask array
    % against the dataArray, this dataArray is the actual scalar
    % array on which the averaging over ellipsoidal area is expected
    % ====================  outputs =============================
    % avg is the average over ellipse area calculated from dataArray
    % totalElem is the total number of pixels in ellipse
    % jvector_x, jvector_y are java Vectors of x,y coordinates
    % of all the interior points of the ellipse
    
    jvector_x = java.util.Vector;
    jvector_y = java.util.Vector;    
    phi = [0,pi/2, pi, 2*pi*3/4.0];
    cosphi = cos(phi);
    sinphi = sin(phi);
    theta = pi* orientation/180;   
    R = [ cos(theta)   sin(theta)
         -sin(theta)   cos(theta)];
    % calculate the four end points of major and minor axes 
    xy = [a*cosphi; b*sinphi];
    % generate points on major axis
    xx = ceil(min(xy(1,:))):floor(max(xy(1,:)));
    totalElem = 0;
    sumTotal = 0;
    % for each point on major axis calculate all possible
    % y locations that fall within ellipse
    for i = xx
        proot = [-sqrt((1-i*i/a/a)*b*b), sqrt((1-i*i/a/a)*b*b)];
        yy = real((min(proot))):real((max(proot)));
        xcopy = ones(1,length(yy)).*i;
        % rotate the ellipse that was generated around (0,0) point 
        % to match the oreintation given and translate the ellipse  to the
        % centroid location
        rtxy = R*[xcopy;yy];
        xcor = (rtxy(1,:)+cx);
        ycor = (rtxy(2,:)+cy);
          % establish bound check
        xout1 = find(xcor < 1);
        xout2 = find(xcor > size(bw,2));
        xcor([xout1, xout2]) = [];
        ycor([xout1, xout2]) = [];
        yout1 = find(ycor < 1);
        yout2 = find(ycor > size(bw,1));
        xcor([yout1, yout2]) = [];
        ycor([yout1, yout2]) = [];
        xcor = round(xcor);
        ycor = round(ycor);
        % delete all the points that fall on the background
        % of the binary mask image
        if ~isempty(xcor)
            ii = length(xcor);
            jj = 0;
            while ii  % pseudo condition
                jj = jj + 1;
                if(jj > ii), break; end
                if(bw(ycor(jj),xcor(jj)) < 1)
                    xcor(jj) = [];
                    ycor(jj) = [];
                    ii = ii - 1;
                end            
            end
        end
        if ~isempty(xcor)
            totalElem = totalElem + numel(xcor);
            for ii = 1:length(xcor)
                % xcor has been calculated so far in column direction
                % ycor has been calcuated in row direction
                % remeber graph x direction is actually column direction
                % and graph y direction is row direction 
                sumTotal = sumTotal + dataArray(ycor(ii),xcor(ii));
                jvector_x.addElement(xcor(ii));
                jvector_y.addElement(ycor(ii));
            end
        end
    end
    avg = sumTotal/totalElem;
    totalElem = round(pi*a*b);
end

Calling A Fortran Subroutine from Python

Calling fortran subroutine from python is really easy using the f2py tool that does the compiling and wrapping of fortran code and ultimately produces a python callable module. The process is really painless and fun! Now the question comes why you would want to call fortran from python.  One of the important reasons could be plotting ! Yes, if you had ever wanted to do all the post processing stuffs apart from the main simulation in a language that provides easy-to-use plotting libraries you would have possibly ruled out all basic languages like fortran/cpp/c/java. Most of the people I know either use matlab or python+numpy for post processing for obvious reasons. Well, life does not go smooth and we can not get all the goodies all the time. There might come a time when you may feel frustrated by the speed and memory management of python or matlab. In such situations binding fortran with python may provide you with a relief. In my case, I was trying to calculate correlation between different x-y planes of a large 3D array having dimensions or ranks 2048x2048x64. The data file takes approximately of space. I was trying to calculate correlation using Matlab's xcorr2 and it was excruciatingly slow. It is not like Matlab is slow. I picked the wrong function for my purpose. For large arrays xcorr2 or scipy.signal's 'correlate2d' is not the right choice for correlation calculations and I am going to show why. Before that, lets see how to call fortran from python. The key tool for this purpose is f2py. I use  Debian/Ubuntu linux systems and it was installed from the repository. To get desired output and avoid surprise you just  need to declare variables in fortran subroutines with the 'intent' attribute. From 'intent(in)' or 'intent(out)' attribute f2py will figure out how to arrange the dummy variables in call signature. f2py removes all variables with 'intent(out)' attribute. The resulting module from f2py will return the original 'intent(out)'  variables in a tuple when called from python. By the way I have not used any deferred shape arrays in fortran subroutine dummy arguments.  The fortran file that I compiled using f2py is "multiLevelCrossCorrFFTM.F90" and created a module named "mxcorr". "multiLevelCrossCorrFFTM.F90" uses the  "fftw3" library for calculating correlation. libfftw3.so will be dynamically linked to f2py generated module without any issue. However I had hard time forcing the compiler to search for 'fftw.h' header file. After searching the web I realized that if the subroutine is saved with ".f90" extension the preprocessor does not run correctly to include the header but saving fortran files with " .F90"  extension solves the problem! The produced mxcorr module by f2py is  imported in "forFromPy.py" script and used. Full codes are given and the procedure should be apparent from the example.

A very good tutorial on f2py can be found here : http://www.ucs.cam.ac.uk/docs/course-notes/unix-courses/pythonfortran/files/f2py.pdf

Note : It is important to check the call signature of the functions produced by f2py. Issuing the command  help(mxcorr) in python prompt  produces the output:

NAME
mxcorr
FILE
/home/ahsan/MEGA/javaProgramsPostProcessing/NeutralChannel/turbStruc3D_corr3D/testJava/mxcorr.so
DESCRIPTION
This module 'mxcorr' is auto-generated with f2py (version:2).
Functions:
status,res = multilevelcrosscorrfftm(arr3d,template,nx=shape(arr3d,0),ny=shape(arr3d,1),nz=shape(arr3d,2),snx=shape(template,0),sny=shape(template,1))

compare this with the original fortran call signature

subroutine multilevelcrosscorrfftm(status, res, arr3d, nx, ny, nz, template, snx, sny)
    implicit none    
    include 'fftw3.f'
    integer,intent(in) ::  nx, ny, nz, snx, sny 
    integer, intent(out) :: status
    real(kind=8),dimension(nx,ny,nz), intent(in) :: arr3d
    real(kind=8),dimension(nx,ny,nz), intent(out) :: res

It is obvious that f2py has modified the call signature. Not only it has removed the "intent(out)" dummy arguments from the call signature but also it has rearranged the position of input arguments. It has moved 'template' argument from 7 th position to 2nd (next to arr3d)! Also note 'intent(out)' variables (status, res) are being returned in a tuple in the order they appeared in original fortran subroutine.

file : forFromPy.py

[code language="python"]import numpy as np
from scipy import signal
import mxcorr
import time
import matplotlib.pyplot as plt

nx = 200; ny = 200;
imgf = np.random.random((nx,ny,1))
temf = imgf[:,:,0]
stime = time.time()
status, corrf = mxcorr.multilevelcrosscorrfftm(imgf,temf,nx,ny,1,nx,ny)
etime = time.time()
print "Elapsed time, fortran call :", etime-stime," sec";

corrp = np.zeros((nx,ny));
stime = time.time()
corrp[:,:] = signal.correlate2d(imgf[:,:,0], temf, mode='same')
etime = time.time() 
print "Elapsed time, python call :", etime-stime, " sec";
corrp = corrp/(sum(sum(temf**2))*sum(sum(imgf**2)))**0.50

plt.subplot(2, 2, 1)
plt.imshow(corrf[:,:,0], cmap='RdBu', interpolation='nearest')
plt.colorbar()
plt.subplot(2, 2, 2)
plt.imshow(corrp[:,:], cmap='RdBu', interpolation='nearest')
plt.title('Correlation results')
plt.colorbar()
plt.show()


[/code]

Elapsed time, fortran call :  0.053929805755  sec

Elapsed time, python call :  5.06791615486  sec

Speed Gain : 101 times !!

Well this is one of those extreme cases! AS has been pointed out, If it is needed to calculate correlation between very large arrays, it will not be wise to use the regular template matching routines like Matlab's 'xcorr2' or in this case 'correlate2d' from scipy.signal. I don not know the exact algorithm that are being used in these functions but I have a feeling that the correlation/convolution is being calculated in this cases using the trivial term by term multiplications. The trivial algorithm "slide the template image over all elements of the target image and multiply term by term and take the sum" will have a computational cost of approximately (nx * ny * nz * (nx * ny+ nx * ny)) whereas, calculating the correlation/convolution using optimized fft libraries like fftw3 will incur a total computational cost of approximately (nz * (nx * ny * log(nx * ny))). So the speed up in the execution time is not surprising at all. A very similar function like "multilevelcrosscorrfftm" can be written in python using pythons fft routines and I guess it will be as fast as the fortran routine. The point is, it is important to realize early what algorithm we are using for our calculations. If our data file is small, say, couple of megabytes we will not notice the performance difference. As soon as data gets large and large, picking the right algorithm becomes a paramount task otherwise, one may find himself in utter frustration like me ! I particularly choose fortran in this case because other calculations involved in my project will make use of extensive loop operations, at least 3 level deep for loops and data set is also large. Both Matlab and Python have got bad reputation regarding loop executions. I am just being cautious from the beginning and writing routines using fortran.

file :  multiLevelCrossCorrFFTM.F90 

[code language="cpp"]
! author : ahsan, cfd.ahsan@outlook.com, Aug 2015
subroutine multilevelcrosscorrfftm(status, res, arr3d, nx, ny, nz, template, snx, sny)
!    this subroutine calcuates the correlation between the template image (template) 
!    and target image (arr3d). the template image is a 2d array and can be smaller than 
!    than any of the x-y plane of the target image array (arr3d)
!    the target image is assumed to be a 3d array and each x-y plane is treated
!    as a separate image. the resulting correlation array has the size of the target
!    array. nx, ny, nz are the number of row (assued x-direction), 
!    column (assumed y-direction), and page(assumed z-direction) of the arr3d array.
!`    snx and sny are the number of row and column of the template image 
    implicit none    
    include 'fftw3.f'
    integer, intent(in)::  nx, ny, nz, snx, sny 
    real(kind=8),dimension(nx,ny,nz), intent(in) :: arr3d
    real(kind=8),dimension(snx,sny), intent(in) :: template
    integer, intent(out) :: status
    real(kind=8),dimension(nx,ny,nz), intent(out) :: res
    real(kind=8), dimension(:,:), allocatable  :: planeData, templateData,xx,yy
    real(kind=8), dimension(:), allocatable  :: onerow, onecol
    complex(kind=8), dimension(:,:), allocatable :: templatefft, planeftt, corr
    integer :: i,j,levelRow, levelCol, templateRow, templateCol, m ,n,nrow,ncol
    integer(kind=8) :: planfd, planbd, planft, planbt    
    
    status = 0 ! status 0 mean success    
! ! add an additional row or column at the end to make data and template size odd
    if(mod(nx,2) == 0 .and. mod(ny,2)==0) then
    levelRow = nx+1; levelCol = ny+1    
    else if(mod(nx,2) == 0 .and. mod(ny,2)==1) then
     levelRow = nx+1; levelCol= ny;     
    else if(mod(nx,2) == 1 .and. mod(ny,2)==0) then
    levelRow = nx; levelCol= ny+1;
    else if(mod(nx,2) == 1 .and. mod(ny,2)==1) then
    levelRow = nx; levelCol = ny;
    end if
    
    if(mod(snx,2) == 0 .and. mod(sny,2)==0) then 
    templateRow = snx+1; templateCol = sny+1;
    else if(mod(snx,2) == 0 .and. mod(sny,2)==1) then
    templateRow = snx+1; templateCol= sny;
    else if(mod(snx,2) == 1 .and. mod(sny,2)==0) then
    templateRow = snx; templateCol = sny+1;
    else if(mod(snx,2) == 1 .and. mod(sny,2)==1) then
    templateRow= snx; templateCol = sny;
    end if
    
    allocate(planeData(levelRow+templateRow -1, levelCol+templateCol-1), &

      STAT = status)
    allocate(templateData(levelRow+templateRow -1, levelCol+templateCol-1), &

      STAT = status)
    allocate(templatefft(size(templateData,1)/2+1, size(templateData,2)), &

      STAT=status)
    allocate(planeftt(size(planeData,1)/2+1, size(planeData,2)), STAT=status)
    allocate(corr(size(planeData,1)/2+1, size(planeData,2)), STAT=status)
    allocate(onerow(size(planeData,1)), onecol(size(planeData,2)))
    allocate(xx(size(planeData,1), ((size(planeData,2))-1)/2))
    allocate(yy((size(planeData,1)-1)/2,size(planeData,2)))
    nrow = size(templateData,1); ncol = size(templateData,2);
    do i=1,nz
        planeData = 0.0
        templateData = 0.0
        xx = 0.0; yy = 0.0;
        planeData(1:size(arr3d(:,:,i),1), 1:size(arr3d(:,:,i),2)) = arr3d(:,:,i)
        templateData(1:size(template(:,:),1), 1:size(template(:,:),2)) = template    
            
            
        call dfftw_plan_dft_r2c_2d(planfd,size(planeData,1),size(planeData,2), &
              planeData,planeftt,fftw_estimate)    
        call dfftw_execute_dft_r2c(planfd, planeData,planeftt)    
        call dfftw_destroy_plan(planfd)
        
        call dfftw_plan_dft_r2c_2d(planft,size(templateData,1), &

            size(templateData,2), templateData,templatefft,fftw_estimate)    

        call dfftw_execute_dft_r2c(planft, templateData,templatefft)    
        call dfftw_destroy_plan(planft)
        
        corr = planeftt * conjg(templatefft)
        call dfftw_plan_dft_c2r_2d(planbd,size(planeData,1), &

          size(planeData,2),corr, planeData, fftw_estimate)

        call dfftw_execute_dft_c2r(planbd,corr, planeData)
        call dfftw_destroy_plan(planbd)
        planeData = planeData / size(planeData)
        ! !     since fftw library does not have a fftshift function
        ! !     manually rearrange the data to keep highest auto 

        ! !     correlation at center
        ! !     resulting correlation matrix output will conform to the matlab xcorr2 result
        onecol = planeData(:,1);
        xx = planeData(:,2:(size(planeData,2)-1)/2+1);
        planeData(:,1:(size(planeData,2)-1)/2)= &
           planeData(:, (size(planeData,2)-1)/2+2:size(planeData,2));
        planeData(:, (size(planeData,2)-1)/2+2:size(planeData,2)) = xx;
        planeData(:,(size(planeData,2)-1)/2+1) = onecol   
     
        onerow = planeData(1,:);
        yy = planeData(2:(size(planeData,2)-1)/2+1, : ) 
        planeData(1:(size(planeData,2)-1)/2, :)= &
            planeData((size(planeData,2)-1)/2+2:size(planeData,2),: ) 
        planeData((size(planeData,2)-1)/2+2:size(planeData,2),:) = yy
        planeData((size(planeData,2)-1)/2+1,:) = onerow
        ! !     normalize the result to get maximum correlation of 1
        planeData = planeData/dsqrt(sum(arr3d(:,:,i)**2)*sum(template**2))
        ! !     return the middle portion of resultant matrix that corresponds 
        ! !     to the size of the input matrix
        res(:,:,i) = planeData((templateRow-1)/2+2:size(planeData,1)- &

            (templateRow-1)/2+1, &
            (templateCol-1)/2+2:size(planeData,2)-(templateCol-1)/2+1)

    end do   

    deallocate(planeData, onecol, onerow,planeftt, templatefft,corr, STAT=status)
    deallocate(templateData, xx, yy, STAT=status)
end subroutine multilevelcrosscorrfftm


[/code]

!! compiling with f2py to make a python callable module
!! f2py -I/usr/include -c -L/usr/lib/x86_64-linux-gnu -lfftw3 -m mxcorr multiLevelCrossCorrFFTM.F90