File Exchange

image thumbnail

Accurate subpixel edge location

version 2.13 (246 KB) by Agustin Trujillo-Pino
Detection of subpixel edges with very high precision in grey level images


Updated 23 Oct 2018

View License

This is the Matlab source code of the sub pixel edge detection method detailed in the paper “Accurate Subpixel Edge Location Based on Partial Area Effect”, published by Elsevier in
The main folder contains the following files:
- subpixelEdges is the detection method. Type 'help subpixelEdges' for more info
- visEdges is the method for displaying detected edges over the image
- subpixelImage creates a high resolution binary image with the detected edges
- subsetEdges extracts a subset of edges that fulfill condition
- example1 is an example of using the method with a synthetic image
- example2 is similar but using a real image
- example3 uses a big image of a printed text captured by a mobile phone camera
- example4 creates a high resolution binary image for the example 1
- example5 is an example of using the method with a high-noise synthetic image
- example6 is similar but using a real angiography
- subpixelStartGUI is a graphical app to test different features of the method
Several slides explaining the method in:

Cite As

Agustin Trujillo-Pino (2019). Accurate subpixel edge location (, MATLAB Central File Exchange. Retrieved .

Comments and Ratings (62)


arnold (view profile)

I converted some of the core functions to mex files a while back, got roughly a 70%+ speed up, so that's really worth it IMHO.
The code remained with my old employer, otherwise I'd post it. My personal license doesn't include coder, sadly.

Sorry Ali. Currently I have only written the Matlab code. Maybe you can try Matlab Coder?

Hi Agustin,

I am very interested in your method of edge detection. I am just wondering do you also have C/C++ implementation of this method?


dimitris dsis

uche Osahor

Hi Zora:
This problem is named edge linking. There are several algorithms to do that, but most of them are oriented at a pixel level. But in this case, you have subpixel infomation, so you could do something more accurate. For example, searching for neighbour edge pixels that share similar orientation can be joined.

Another strategy is that, if you assume that the contour line could be approximated by a known curve, like a parabola or an ellipse, you can use a numerical algorithm like polyfit to find the best curve that fits all your subpixel coordinates.

zora zhang

Hi, Agustion:
Thanks for the code, it's robust and correct! I wanner to get a edge chain in order .Your code gives the sub-pixel position of edge points but they aren't linked in order. How could I do?

zhifan zhang

huahua niu

Jonghoon Ryu

Hammad Faizi


Yes Efstathios. You're right. I will add this sentence in the next revision. Thanks

Hi Agustin,

Thanks for the code. I think in subpixelImage.m there should be a prediction like that radii(radii==-inf) = -1e6;
You have one for the Inf but if -inf appeared there would be an error.

Hi Prasangsha: for each edge pixel, the normalized vector <nx,ny> indicates the normal direction to the edge. And the absolute difference between outer and inner intensities, abs(i1-i0) could be considered as the magnitude of the gradient vector in that point.

How to find the magnitude of the vector?


Andrew (view profile)


Robust, fast and accurate !
Thanks a lot !



arnold (view profile)

Dear Augustin. There is nothing wrong with your function. Coder only expects variables to be declared a certain way.
When I find the time and deem it necessary to speed things up I'll have a closer look. If I'm successful I will email you the result. Again, the goal would only be to speed things up more. You could then attach the c-functions with your code and people who want it to be faster just compile them within matlab, it's no big deal.
The speedup I currently get is around a factor of 2.5 but I had to redefine the edge class with a normal structure and then re-cast this struct as an edgepixel-class in subpixel.m. Not elegant, but it works for now. If I figure out how to do it properly, I'll let you know :)

Again, I find your contribution and paper to be very useful!

Sorry arnold, but I don't have experience with coder. And I couldn't understand the error about undefined 'posision' inside EdgePixel class. It's defined in the EdgePixel.m file.


arnold (view profile)

All in all the function is not slow by any standard but looking at the runtimes of the functions I was certain they could get faster. Having a closer look, I see no immediate approach to use the GPU for further speedup.
However, generating c-code (mex file) using coder should increase performance significantly. I'm no expert but I've tried to pack the finaldetectorIter functions into c files.
Just as an example, packing just lines 167-169 or lines 351-353 in finalDetectorIter1 (rimv2 = ...) as mex speeds up these two assignments by more than 400% which is surprising since they get called A LOT. It would be best to pack the entire loop or function, but coder can't seem to handle the assignments into the edge class, at least I haven't figured that one out yet. Coder gives the error:

Property 'position' is undefined'... and that for all "ep." assignments.


arnold (view profile)

Hi Augustin,
thank you very much for the quick fix. Your work is truly appreciated and if we can use this in work to be published we will cite your paper properly!

Just a hint to anyone having to call the function numerous times. I can recommend getting rid of the addpath in line 90, this line is very often the one to cause the biggest overall compute time.
I also found that it makes way more sense to run the code less often on larger arrays than actually splitting up arrays to smaller blocks - which often is faster with other functions.
I don't have a CUDA machine here but I will try and see if using GPUARRAY can speed things up some more.

Hi Arnold: bug fixed. Download and try again


arnold (view profile)

Hi Augustin,

I've encountered an error occuring only with few out of thousands of images:
Error using .*
Matrix dimensions must agree.

Error in finalDetectorIter1 (line 391)
ep.curv = [ep.curv; 2*c.*n./((1+b.^2).^1.5)];

Error in subpixelEdges (line 94)
edges = finalDetectorIter1(image, threshold, order);

I don't understand where it's coming from but changing the threshold sometimes works.
I've sent you an email with a sample image as *.mat.

Hi Bruno: you're right. Thanks. It was a mistake. RI is the name of the restored image to display


Bruno (view profile)

example6.m, line 16:
maybe you mean imshow(RI/255...
instead of


arnold (view profile)

Hi Agustin,
thanks for adding the subsetEdges method, I'll be using it pretty much every time I use this function.

Hi Arnold: a new method subsetEdges has been included to extract a subset of edges that fulfill condition. You can try now


arnold (view profile)

I knew about the 'condition' functionality of the visEdges function and it's useful for displaying, yes.

However, if one wants to do more than view, having to copy the variable one element at a time seems rather cumbersome to do all the time. An ordinary structure would behave more like expected. I'm no expert in class definition though so I can't comment on how easy/hard it would be to make this more straight-forward.

thanks for your input!

Hi Arnold. To do that you can try the following:
condition = edges.x < 10;
visEdges(edges, 'condition', condition);
edgesC = EdgePixel;
edgesC.x = edges.x(condition);
edgesC.y = edges.y(condition);


arnold (view profile)

Hi Augustin,
how can one easily access the EdgePixel class? Say I'd want to filter it or create a subset via indexing.

edges(edges.x<10) = []

All I get is: " Index exceeds matrix dimensions."

Hi Arnold: could you send me a sample image? Maybe another approaches could be tried.


arnold (view profile)

thank you Agustin. Neither of your proposals offer subpixel accuracy though.
I'm currently working on a detection function of laser marks/fiducials (bright crosses on dark ground) and I would like to have a robust subpixel accuracy included, that's how I stumbled upon this one :)

Hi arnold: actually this algorithm only detects edges. It’s not a ridge detector. So two different edges will be detected in the presence of thick lines. In your case, centerline detection or skeletonization methods would be more appropiate. You can try this:


arnold (view profile)

Hi Augustin,

I'm trying to find bright lines on a dark background. Your algorithm works fine and finds the edges, but I'm more interested in the 'ridge' (the center of each line).
Now I find one edge on each side of the lines, 2px apart. Do you see any obvious path to extract this from your edge vectors? Basically that would be a 'ridge detector'... with gradients on both sides. In your example image that would be the center of the ring, instead of the inner and outer edge.

Hi Gabby: what you want to do is called 'edge linking'. There are several algorithms to do that. There is also several Matlab functions that solve this problem (like edgelink), but is at a pixel level (not subpixel).
But the theory inside an edge linking algorithm is very easy: for every subpixel position (you also has an estimation of the orientation), you have to look for neighboor pixels in that direction (one pixel at each side), and join the subpixel position with the positions of the neighbours. And you've got it.

Jannie Sy

Thanks Agustin. I want to connect the subpixel resuts such that they produce a continuous line that a typical edge detector like Canny does? I tried the canny and sobel operators but I am not getting enough detail. When I saw this, it produced more accurate results when it comes to edge detection. However, the results are subpixel points only. How can I make them continuous?

Hi Gabby: depends on what kind of line you expect to obtain:
- if straigth line, you can use polyfit with degree one,
- if circle line, you can use circfit
- if looking for a polyline, you can adjust all the neighbour points that be colinear to small straight segments
- for curves, you can use polyfit with higher degrees,

Jannie Sy

Hi, thank you for this. I have a question, do you any way I can connect the points (the subpixel results) in order to make a line? I tried roipoly but I'm having a difficulty in it.

Hi Siew Hon:

Effectively a bug was found and fixed. Download the new version and try again please.

Siew Hon Teay


Thanks for sharing the code. I try to implement the algorithm to detect edge for a sub-image of a bearing image. However, this is the error that I got:

Subscript indices must either be real positive integers or logicals.

Error in subpixelImage (line 40)
bw((round(x)-1)*rows+round(y)) = true;

Could you please guide on fixing this issue? Thanks!

Enclosed please find the link for the image:

Hi Xin Gao: could you share your image so that I can reproduce this error?

Xin Gao

I use an 480*720 image to perform edge detection (R2015a), %% subpixel detection
threshold = 16;
edges = subpixelEdges(image, threshold);
%% show edges
visEdges(edges, 'showNormals', false);
binaryImage = subpixelImage(image, size(image), 1);
% figure; imshow(binaryImage, 'InitialMagnification', 'fit');
however, it reports the following error:
Attempt to reference field of non-structure array.
Error in subpixelImage (line 14)
radii = 1 ./ edges.curv; % Problem is located here!
Error in example2 (line 18)
binaryImage = subpixelImage(image, size(image), 1);

Anybody expert can fix this error?
Thank you very much!

That is really kind of you. Thanks!
It's a great code.
I hope I can run this on an FPGA, since that is where it has to perform ultimately.

Hi Meghana:

This method just detects all the pixels belonging to edges, and computes all its subpixel features. But it does not produce any continous line. You would need a second algorithm to do that, starting from these subpixel results.

Anyway, I have worked on your problem and estimated a result for the vertical line you are finding, inside the same web page where you made your question:


Thanks for sharing this code.

I have a query. For my application, the edges will always be vertical (almost). An example image is posted here:

This is how I am performing edge detection using your functions:

image = imread('img1.bmp');
%% subpixel detection
disp('Computing subpixel edges...');
threshold = 30;
edges = subpixelEdges(image, threshold);
%% show edges
visEdges(edges, 'showNormals', false);
bw = subpixelImage(edges, size(image), 10);
figure; imshow(bw, 'InitialMagnification', 'fit');

The output of this seems to be dis-continuous.
Is there any way by which I can get one single continuous line? (maybe by adjusting different parameters) Is this within the scope of this function?


Compan S2

Very fast and accurate



Xen (view profile)

Just to add, this method is also ~8 times faster.


Xen (view profile)

Works great Agustin. I compared your method with Canny edge detection on a resized image (by bicubic interpolation) and I get pretty similar results, with the added advantage of your method having just one parameter to optimize.

Hi Xen:
Yes, it would be possible. Let's suppose a new binary 5x image, where every original pixel corresponds to a grid of 5x5 small pixels in the binary image. For every pixel edge, as you have all the edge parameters, a second order curve could be plotted over the binary image small pixels.
I have added this new method, named subpixelImage(). Just try example4 to test it.


Xen (view profile)

Great work Agustin.
Is it possible to save the edges in a simple binary image? And I mean, of course, in a higher resolution than the original input image (e.g. 5x or 10x larger), otherwise it wouldn't be much useful for sub-pixel analysis!

Hi Tolga:

In our experiments, the averaging kernel is more effective to remove noise than gaussian. Now the source code is only defined for this mask. But we can add the option to use a specific smooth mask. And for the rest of your suggestions, I agree with you that using conv2 and gradient methods will produce a more clear code. Both will be added in the next version.

As for your second question, we are not interested in the specific partial derivative values as in the differences of intensity values inside the same row / column. That’s the reason why we use only Gx and Gy

Tolga Birdal


Thank you for the submission. To clarify the implementation, I have found the following to be useful:

1. Instead of an averaging kernel, one could use a Gaussian one:

g = fspecial('gaussian');
G = conv2(F,g,'same');

2. Your gradient computation is the same as:
[Gx, Gy] = gradient(G);
grad = sqrt(Gx.^2+Gy.^2);


Tolga Birdal

I guess MATLAB's gradient function would be a replacement for the computation of the partial derivatives. Yet, in many cases Sobel is a more robust approximation. Is there any reson that you are not reverting to those?



Minor changes

Minor updates
Fixed references to external test images

Update link with slides explaining the algorithm

Bug solved when non valid edges detected

New method subsetEdges to extract a subset of edges that fulfill condition

New option to allow several iterations of the method (smoothing + detection + synthetic image creation) for better subpixel edge detection in high-noise images.

Added link to slides explaining the algorithm

Bug fixed when invalid values are obtained

New example4 added

Mistakenly marked as toolbox

New method subpixelImage to compute a high resolution binary image with the detected edges

New method subpixelImage to compute a high resolution binary image with the detected edges

New description

Update directories belonging on the MATLAB Search Path

MATLAB Release Compatibility
Created with R2013a
Compatible with any release
Platform Compatibility
Windows macOS Linux

Discover Live Editor

Create scripts with code, output, and formatted text in a single executable document.

Learn About Live Editor