Thread Subject: Need algorithm help: molecule length from an image

Hello,

My name is Adrian Randall and I work in an Electron
Microscopy lab. I am trying to determine the lengths of
molecules from images using MATLAB. I have just finished a
semester of MATLAB, so I'm very familiar with the program. I
need help coming up with an algorithm for such a program.
Does anyone have any experience with image analysis that
could give me a hand?

Here is a link to what an electron micrograph looks like:
http://www.unc.edu/~arandall/image.png

Thank you in advance for any suggestions!

-Adrian Randall (ajr [at] med.unc.edu)

In article <g023qb$j0m$1@fred.mathworks.com>, Adrian <ajr@med.unc.edu> wrote:

>My name is Adrian Randall and I work in an Electron
>Microscopy lab. I am trying to determine the lengths of
>molecules from images using MATLAB.

I had a glance at the image, and the first question in my mind
was whether you want the length as folded (e.g., the length
of an appropriate oriented bounded box), or if you want the
length as if were unfolded (length of the perimeter) ? If as if
unfolded, then how should the loop at one end be treated, as
a circle whose diameter is to be taken, or the circumference
(af if the loop were stretched out to its maximum length) ?
--
  "Man's life is but a jest,
   A dream, a shadow, bubble, air, a vapor at the best."
                                         -- George Walter Thornbury

On May 9, 2:06=A0pm, "Adrian " <a...@med.unc.edu> wrote:
> Hello,
>
> My name is Adrian Randall and I work in an Electron
> Microscopy lab. I am trying to determine the lengths of
> molecules from images using MATLAB. I have just finished a
> semester of MATLAB, so I'm very familiar with the program. I
> need help coming up with an algorithm for such a program.
> Does anyone have any experience with image analysis that
> could give me a hand?
>
> Here is a link to what an electron micrograph looks like:http://www.unc.ed=
u/~arandall/image.png
>
> Thank you in advance for any suggestions!
>
> -Adrian Randall (ajr [at] med.unc.edu)

Adrian:
Thanks for posting the image in advance (that's one better than most
people). I do have nearly 30 years of image analysis experience and
might be able to help you but you didn't say what you were looking for
and I don't recognize what is what in the image. Is this an atomic
force microscopy image (the only one that I know of that allows you to
see individual atoms)? Can you have the system imprint a scale bar
onto the image? You need to do a spatial calibration on the image so
to determine the length of a molecule you first need some sort of
scale bar or other way of knowing how many pixels equal how many real
world unit. If that molecule is that long squiggle then there are
many ways of finding it. You could even just manually trace it. But
then when you have the number of pixels long it is, you'll want to
convert than into nanometers or something. This is where you'll need
the spatial calibration.

To find the squiggle, you could try something like a big blurring
(convolution) filter and subtract the original image, threshold, then
clean up a bit (maybe by doing a morphological closing operation to
make sure the squiggle is not in multiple pieces or a smaller median
filter to throw out small particles), then bwlabel() and regionprops()
and toss out all blobs except the blob with the largest area. Then
look at its perimeter property. Then multiply by your spatial
calibration factor. That's what I'd try first and then I'd modify the
algorithm depending on the quality of results I achieved.
Regards,
ImageAnalyst

"Adrian " <ajr@med.unc.edu> wrote in message <g023qb$j0m
$1@fred.mathworks.com>...
> Hello,
>
> My name is Adrian Randall and I work in an Electron
> Microscopy lab. I am trying to determine the lengths of
> molecules from images using MATLAB. I have just finished a
> semester of MATLAB, so I'm very familiar with the program. I
> need help coming up with an algorithm for such a program.
> Does anyone have any experience with image analysis that
> could give me a hand?
>
> Here is a link to what an electron micrograph looks like:
> http://www.unc.edu/~arandall/image.png
>
> Thank you in advance for any suggestions!
>
> -Adrian Randall (ajr [at] med.unc.edu)
------------
  Adrian, there is no denying that your problem is a difficult one and involves
several stages. In my opinion the most difficult stage comes after you have
successfully processed an image to the point where it is a single connected
black and white image without any false holes within it, giving a good
silhouette of the laid-out molecule. What I describe here is by no means a
solution but merely an insight into what kinds of problems need to be solved.

  Think how human beings would proceed if faced with this problem. They
would probably start sketching little line segments or arcs which were more
or less tangent to the shape at various points, and then they would begin
drawing a kind of envelope along these lines or curves. This would hopefully
produce a single curve running along the centered-length of the shape.
Finally they would probably whip out their handy little arc length tool with the
calibrated wheel at one end and roll it along the total length of the curve to
get its length.

  This should give you an idea of the formidable tasks that await you in this
project if you want to develop a fully automated and reliable answer to your
question. You will need to find best-fitting line segments or arcs to various
regions of the shape. You need to be able to connect these into a single
coherent envelope curve along the length of the shape. Finally you need to
determine the length of this envelope.

  Note to ImageAnalyst: In my opinion the length of the shape's perimeter
(presumably divided by 2) is not a reliable measure of what Adrian is seeking
as the over-all length of a molecule. Its outline is likely to be a somewhat
jagged affair as is characteristic of long-chain molecules with differing kinds
of structures and would tend to give too large and inconsistent a value. I
think his problem is more difficult than that.

Roger Stafford

In article <g07jfi$afk$1@fred.mathworks.com>,
Roger Stafford <ellieandrogerxyzzy@mindspring.com.invalid> wrote:

> Note to ImageAnalyst: In my opinion the length of the shape's perimeter
>(presumably divided by 2) is not a reliable measure of what Adrian is seeking
>as the over-all length of a molecule. Its outline is likely to be a somewhat
>jagged affair as is characteristic of long-chain molecules with differing kinds
>of structures and would tend to give too large and inconsistent a value.

Hmmm, that triggers the throught that the molecule may have significant
3 dimensional components; unless we have at least two views of it
in different planes, we would not be able to tell whether the molecule
just happens to be shorter in the direction of view, but long into
or out of the direction of view.
--
  "The quirks and arbitrariness we observe force us to the
  conclusion that ours is not the only universe." -- Walter Kistler

On May 11, 7:12=A0pm, rober...@ibd.nrc-cnrc.gc.ca (Walter Roberson)
wrote:
> In article <g07jfi$af...@fred.mathworks.com>,
>
> Roger Stafford <ellieandrogerxy...@mindspring.com.invalid> wrote:
> > =A0Note to ImageAnalyst: In my opinion the length of the shape's perimet=
er
> >(presumably divided by 2) is not a reliable measure of what Adrian is see=
king
> >as the over-all length of a molecule. =A0Its outline is likely to be a so=
mewhat
> >jagged affair as is characteristic of long-chain molecules with differing=
 kinds
> >of structures and would tend to give too large and inconsistent a value.
>
> Hmmm, that triggers the throught that the molecule may have significant
> 3 dimensional components; unless we have at least two views of it
> in different planes, we would not be able to tell whether the molecule
> just happens to be shorter in the direction of view, but long into
> or out of the direction of view.
> --
> =A0 "The quirks and arbitrariness we observe force us to the
> =A0 conclusion that ours is not the only universe." -- Walter Kistler

----------------------------------------------------------------------------=
----------------
Well Walter, I believe you're exactly right. I imagine that, assuming
the molecule is the string-like object in the scene, that if you took
a picture of it at one time, at another time the molecule could have
squiggled to a completely different shape, and since it's probably
floating around in 3D in some liquid, its projected shape (onto a 2D
image) and length could be completely different from one image to the
next. This change is called "repeatability" and is how this one
molecule specimen changes from one snap of that one molecule to the
next. There's also something called "reproducibility" which is a
different specimen of the same species of molecule but imaged as a
separate sample on perhaps a different microscope. So because of all
this, it's possible that the molecule could have different "lengths"
depending on when you snapped the picture because you're seeing just a
2D projection of it.

And another thing . . . when you take a look at the image he posted,
let's say the "head" of the molecule is that wad of curled up stuff
plus the loop, and the tail is the other end. Now exactly where does
the tail end? You can see that it's not so easy to answer. Maybe it
curls around down below and heads back to the right. (It's not as
obvious now - I think he replaced the with a smaller version of what
he used to have). But the tail seems to sort of fade away and it's
not really clear where it ends. It is sort of a judgement call. And
what about the head? Is there a bunch of length there that is all
tangled up into a ball so that you can't really know the "true"
length?

And Adrian didn't say what kind of precision he needs. Maybe a fully
automated method with sub-pixel precision is needed, but maybe just
hand-tracing the line is good enough especially if the length changes
a lot from one snapshot to the next.

Because of all this, it probably doesn't make sense to go overboard
with fancy algorithms to get the length of the molecule in that one
single, specific picture with sub-pixel precision when the "length" of
the molecule may vary by tens or hundreds of pixels from one picture
to the next.

Hey Adrian, are you ever going to check back here for answers to your
question, and answer some of the issues brought up?
Regards,
ImageAnalyst

ImageAnalyst <imageanalyst@mailinator.com> wrote in message
<7d6b1108-5156-49cb-84bf-
423bb11b437a@k13g2000hse.googlegroups.com>...
> On May 11, 7:12=A0pm, rober...@ibd.nrc-cnrc.gc.ca (Walter
Roberson)
> wrote:
> > In article <g07jfi$af...@fred.mathworks.com>,
> >
> > Roger Stafford
<ellieandrogerxy...@mindspring.com.invalid> wrote:
> > > =A0Note to ImageAnalyst: In my opinion the length of
the shape's perimet=
> er
> > >(presumably divided by 2) is not a reliable measure of
what Adrian is see=
> king
> > >as the over-all length of a molecule. =A0Its outline
is likely to be a so=
> mewhat
> > >jagged affair as is characteristic of long-chain
molecules with differing=
> kinds
> > >of structures and would tend to give too large and
inconsistent a value.
> >
> > Hmmm, that triggers the throught that the molecule may
have significant
> > 3 dimensional components; unless we have at least two
views of it
> > in different planes, we would not be able to tell
whether the molecule
> > just happens to be shorter in the direction of view,
but long into
> > or out of the direction of view.
> > --
> > =A0 "The quirks and arbitrariness we observe force us
to the
> > =A0 conclusion that ours is not the only universe." --
Walter Kistler
>
> ----------------------------------------------------------
------------------=
> ----------------
> Well Walter, I believe you're exactly right. I imagine
that, assuming
> the molecule is the string-like object in the scene, that
if you took
> a picture of it at one time, at another time the molecule
could have
> squiggled to a completely different shape, and since it's
probably
> floating around in 3D in some liquid, its projected shape
(onto a 2D
> image) and length could be completely different from one
image to the
> next. This change is called "repeatability" and is how
this one
> molecule specimen changes from one snap of that one
molecule to the
> next. There's also something called "reproducibility"
which is a
> different specimen of the same species of molecule but
imaged as a
> separate sample on perhaps a different microscope. So
because of all
> this, it's possible that the molecule could have
different "lengths"
> depending on when you snapped the picture because you're
seeing just a
> 2D projection of it.
>
> And another thing . . . when you take a look at the image
he posted,
> let's say the "head" of the molecule is that wad of
curled up stuff
> plus the loop, and the tail is the other end. Now
exactly where does
> the tail end? You can see that it's not so easy to
answer. Maybe it
> curls around down below and heads back to the right.
(It's not as
> obvious now - I think he replaced the with a smaller
version of what
> he used to have). But the tail seems to sort of fade
away and it's
> not really clear where it ends. It is sort of a
judgement call. And
> what about the head? Is there a bunch of length there
that is all
> tangled up into a ball so that you can't really know
the "true"
> length?
>
> And Adrian didn't say what kind of precision he needs.
Maybe a fully
> automated method with sub-pixel precision is needed, but
maybe just
> hand-tracing the line is good enough especially if the
length changes
> a lot from one snapshot to the next.
>
> Because of all this, it probably doesn't make sense to go
overboard
> with fancy algorithms to get the length of the molecule
in that one
> single, specific picture with sub-pixel precision when
the "length" of
> the molecule may vary by tens or hundreds of pixels from
one picture
> to the next.
>
> Hey Adrian, are you ever going to check back here for
answers to your
> question, and answer some of the issues brought up?
> Regards,
> ImageAnalyst

I would be interested in Adrian's comments too but I've
enjoyed following the thread and thinking about the
different ways of approaching the problem. There have been
some good ideas.

I am wary of answering questions that Adrian really needs
to answer but I think the image shows a molecule that has
been deposited on a surface and then stained with something
to increase the contrast for electron microscopy. This
means that the molecule is deposited 'flat' on a surface
and doesn't cross over itself. I agree that the problem is
harder if this isn't the case. For what it's worth I don't
think that you can do electron microscopy on molecules 'in
solution'.

I've seen people talk about similar images in seminars I've
attended and I'd be surprised if Adrian needs to 'unravel'
the ball or requires 'sub-pixel' methods. I think
the 'perimeter' approach is probably the most appropriate
idea that's been suggested so far for this type of work.

But as ImageAnalyst points out, that's really for Adrian to
decide.

Ken

Thank you all for your insight into the matter. Sorry for
the lack of clarification for the image I linked to: the
squiggly line in the middle is in fact the molecule.
Everything else is background. The image is a transmission
electron micrograph, so it's 2D and fixed.

The way that we do length analysis is by manually drawing
along the molecule using a mouse and some software that came
with the electron microscope. As you can image, this is
somewhat inaccurate and tedious, especially if one is
scoring hundreds of molecules. In addition, for the project
I am working on, length analysis is the centerpiece, so I
want to have data as accurate as possible.

I didn't include a scale bar - my apologies. I wanted to
isolate only one molecule to make the problem easier to
solve, and forgot to include the scale. There is a scale
that converts pixel length to actual distance that I know
how to use, so the only problem is finding the pixel length.

This is a method I have come up with from using your
suggestions.

1) Find the average pixel shade (between 0 and 255)
2) Remove all pixels that are below the threshold (remove
the darker background).
3) Set the remaining pixels equal to 1; all other pixels
equal to 0.
4) Scan the image from left to right and then "walk" along
areas that are labeled as "1". If the program finds a pixel
that is labeled "1", look toward the 8 surrounding pixels to
find if there is another pixel labeled "1". If so, continue
to look for pixels labeled "1" (e.g., keep "walking" along
the molecule), and if possible, continue "walking" in the
same direction as previously (to prevent the "walker" from
going back to where it came from). Then, after walking along
the molecule, label all of the points it walked on as "0" so
the scanner won't walk there again.
5) Output the length measurement as pixel and convert to nm
by looking at the scale.

Does this seem reasonable? Thank you very much for the
discussion..

On May 13, 3:08=A0pm, "Adrian " <a...@med.unc.edu> wrote:
> Thank you all for your insight into the matter. Sorry for
> the lack of clarification for the image I linked to: the
> squiggly line in the middle is in fact the molecule.
> Everything else is background. The image is a transmission
> electron micrograph, so it's 2D and fixed.
>
> The way that we do length analysis is by manually drawing
> along the molecule using a mouse and some software that came
> with the electron microscope. As you can image, this is
> somewhat inaccurate and tedious, especially if one is
> scoring hundreds of molecules. In addition, for the project
> I am working on, length analysis is the centerpiece, so I
> want to have data as accurate as possible.
>
> I didn't include a scale bar - my apologies. I wanted to
> isolate only one molecule to make the problem easier to
> solve, =A0and forgot to include the scale. There is a scale
> that converts pixel length to actual distance that I know
> how to use, so the only problem is finding the pixel length.
>
> This is a method I have come up with from using your
> suggestions.
>
> 1) Find the average pixel shade (between 0 and 255)
> 2) Remove all pixels that are below the threshold (remove
> the darker background).
> 3) Set the remaining pixels equal to 1; all other pixels
> equal to 0.
> 4) Scan the image from left to right and then "walk" along
> areas that are labeled as "1". If the program finds a pixel
> that is labeled "1", look toward the 8 surrounding pixels to
> find if there is another pixel labeled "1". If so, continue
> to look for pixels labeled "1" (e.g., keep "walking" along
> the molecule), and if possible, continue "walking" in the
> same direction as previously (to prevent the "walker" from
> going back to where it came from). Then, after walking along
> the molecule, label all of the points it walked on as "0" so
> the scanner won't walk there again.
> 5) Output the length measurement as pixel and convert to nm
> by looking at the scale.
>
> Does this seem reasonable? Thank you very much for the
> discussion..

--------------------------------------------------------------------------
Adrian:
Glad you're still monitoring and we're not just wasting our time
here. You could do what you described with this algorithm:
1. Binarize the image at some threshold.
2. Skeletonize the image - bwmorph(thresholdedImageArray, 'skel',
Inf). That's what I'm doing but you're essentially doing bwlabel and
then binarizing the labeled image. Skeletonizing has the advantage
that it will ignore parts of the line that are thicker than one pixel
but it will under count big solid balled up areas. Your method will
overcount parts of the line that are thicker than one pixel but will
correctly count big balled up areas.
3. Sum the pixels in the image Sum(sum(skeletonizedImageArray))
The only point I might add is that you might want to count triple and
quad points (where the line touches or crosses itself) twice rather
than avoiding counting those pixels like you said. The molecule would
actually occur twice at those positions. Think about it and you'll
see why.
Regards,
ImageAnalyst

In article <g0couj$np$1@fred.mathworks.com>, Adrian <ajr@med.unc.edu> wrote:

>4) Scan the image from left to right and then "walk" along
>areas that are labeled as "1". If the program finds a pixel
>that is labeled "1", look toward the 8 surrounding pixels to
>find if there is another pixel labeled "1". If so, continue
>to look for pixels labeled "1" (e.g., keep "walking" along
>the molecule), and if possible, continue "walking" in the
>same direction as previously (to prevent the "walker" from
>going back to where it came from). Then, after walking along
>the molecule, label all of the points it walked on as "0" so
>the scanner won't walk there again.
>5) Output the length measurement as pixel and convert to nm
>by looking at the scale.

Suppose you have

1011
0100
0010

And you start at the top-left corner. You walk left to right so you
encounter the top-left 1. Presumably you increment your count to 1.
You examine the surrounding pixels and find
the 1 diagonally down and right. You try to continue left to right, but
there is no 1 there so you take the path that is open by going diagonally
down and right. However, as you were not able to continue in the same
direction as you had been going, you label the points you walked in
that stretch with a 0 so the scanner won't walk there again, leading to

0011
0100
0010

You are at the second pixel, presumably you increment your count to 2.
Now you are headed diagonally down and right. You examine the
surrounding pixels and find the 1 diagonally up and right and
the 1 diagonally down and right. As per your algorithm, you prefer
to keep going in the same direction, so you take the diagonal
down and right. As you -were- able to continue in the same direction,
you are on a run and so do -not- mark the pixel you came from with a 0.

You are at the third pixel, presumably you increment your count to 3.
You look around from that bottom pixel and find the only path open
to be the one up and left. However, as you were not able to continue
in the same direction as you had been going, you label the points
you walked in the run with a 0 won't walk there again, leading to

0011
0000
0000

There is now an ambiguity in your algorithm: you had an available
path up-left a moment ago, but that pixel is now 0'd, and there
are now no pixels left in the neighbourhood that are non-0. Do you
stop at this point with a pixel count of 3? My interpretation of
your algorithm is that you do -not- stop, because your instruction
was to do the search before the 0'ing.

So you are now back at 2nd row 2nd column. Your algorithm didn't
say anything about not counting on backtracks, so presumably you
increment the count to 4. You look around and find the open path
to the up-right. That's a change of direction, so you 0 out the pixels
you had in that run, which makes no change because they are already 0.

You are now at the 1st row 3rd column, trying to head up and right.
Presumably you increment your counter to 5. You look around and
find the open path to the right. That's a change of direction
so you 0 out the pixels you had in that run, which makes the matrix

0001
0000
0000

You are now at the 1st row 4th column, trying to head right.
Presumably you increment your counter to 6. You look around and find
no open paths. Presumably you stop there.

Your recorded molecule length is now 6, but the longest connected
chain was 4 -- (1,1) -> (2,2) -> (1,3) -> (1,4), so arguably the
molecule length is 4 instead. Or, since there were two diagonal
moves and then 1 lateral move, 2*sqrt(2)+1, which is about 3.8.
Add 1 if you are counting from the edges instead of from the centers.


Now, suppose you had started with

101111
010000
001000
000100

then when you are examining (2,2) after coming from (1,1),
because your algorithm always continues in the same direction if it
can, it would continue on the downward diagonal to the edge, find
the upward-left step from the bottom as the only connection to the
bottom pixel, and would 0 out the entire diagonal, leaving you with

001111
000000
000000
000000

with you positioned at (3,3). There are no connections from there
so the algorithm would stop, unable to reach that chain at the top
right corner. The algorithm has thus clearly not found the proper length.


Generally speaking, the algorithm's preference for continuing in
the same direction as possible is going to create problems for you,
as it will be unable to recognize where the molecule changes
direction even just due to a standard convex curve. And as I showed
above, counting pixels visited is not an accurate length measure
if any branches exist. Thirdly, you are not talking into account
that diagonal moves are further than vertical or horizontal moves.
--
  "Prevention is the daughter of intelligence."
                                              -- Sir Walter Raleigh

roberson@ibd.nrc-cnrc.gc.ca (Walter Roberson) wrote in
message <g0crs0$6al$1@canopus.cc.umanitoba.ca>...
> In article <g0couj$np$1@fred.mathworks.com>, Adrian
<ajr@med.unc.edu> wrote:
>
> >4) Scan the image from left to right and then "walk"
along
> >areas that are labeled as "1". If the program finds a
pixel
> >that is labeled "1", look toward the 8 surrounding
pixels to
> >find if there is another pixel labeled "1". If so,
continue
> >to look for pixels labeled "1" (e.g., keep "walking"
along
> >the molecule), and if possible, continue "walking" in the
> >same direction as previously (to prevent the "walker"
from
> >going back to where it came from). Then, after walking
along
> >the molecule, label all of the points it walked on
as "0" so
> >the scanner won't walk there again.
> >5) Output the length measurement as pixel and convert to
nm
> >by looking at the scale.
>
> Suppose you have
>
> 1011
> 0100
> 0010
>
> And you start at the top-left corner. You walk left to
right so you
> encounter the top-left 1. Presumably you increment your
count to 1.
> You examine the surrounding pixels and find
> the 1 diagonally down and right. You try to continue left
to right, but
> there is no 1 there so you take the path that is open by
going diagonally
> down and right. However, as you were not able to continue
in the same
> direction as you had been going, you label the points you
walked in
> that stretch with a 0 so the scanner won't walk there
again, leading to
>
> 0011
> 0100
> 0010
>
> You are at the second pixel, presumably you increment
your count to 2.
> Now you are headed diagonally down and right. You examine
the
> surrounding pixels and find the 1 diagonally up and right
and
> the 1 diagonally down and right. As per your algorithm,
you prefer
> to keep going in the same direction, so you take the
diagonal
> down and right. As you -were- able to continue in the
same direction,
> you are on a run and so do -not- mark the pixel you came
from with a 0.
>
> You are at the third pixel, presumably you increment your
count to 3.
> You look around from that bottom pixel and find the only
path open
> to be the one up and left. However, as you were not able
to continue
> in the same direction as you had been going, you label
the points
> you walked in the run with a 0 won't walk there again,
leading to
>
> 0011
> 0000
> 0000
>
> There is now an ambiguity in your algorithm: you had an
available
> path up-left a moment ago, but that pixel is now 0'd, and
there
> are now no pixels left in the neighbourhood that are non-
0. Do you
> stop at this point with a pixel count of 3? My
interpretation of
> your algorithm is that you do -not- stop, because your
instruction
> was to do the search before the 0'ing.
>
> So you are now back at 2nd row 2nd column. Your algorithm
didn't
> say anything about not counting on backtracks, so
presumably you
> increment the count to 4. You look around and find the
open path
> to the up-right. That's a change of direction, so you 0
out the pixels
> you had in that run, which makes no change because they
are already 0.
>
> You are now at the 1st row 3rd column, trying to head up
and right.
> Presumably you increment your counter to 5. You look
around and
> find the open path to the right. That's a change of
direction
> so you 0 out the pixels you had in that run, which makes
the matrix
>
> 0001
> 0000
> 0000
>
> You are now at the 1st row 4th column, trying to head
right.
> Presumably you increment your counter to 6. You look
around and find
> no open paths. Presumably you stop there.
>
> Your recorded molecule length is now 6, but the longest
connected
> chain was 4 -- (1,1) -> (2,2) -> (1,3) -> (1,4), so
arguably the
> molecule length is 4 instead. Or, since there were two
diagonal
> moves and then 1 lateral move, 2*sqrt(2)+1, which is
about 3.8.
> Add 1 if you are counting from the edges instead of from
the centers.
>
>
> Now, suppose you had started with
>
> 101111
> 010000
> 001000
> 000100
>
> then when you are examining (2,2) after coming from (1,1),
> because your algorithm always continues in the same
direction if it
> can, it would continue on the downward diagonal to the
edge, find
> the upward-left step from the bottom as the only
connection to the
> bottom pixel, and would 0 out the entire diagonal,
leaving you with
>
> 001111
> 000000
> 000000
> 000000
>
> with you positioned at (3,3). There are no connections
from there
> so the algorithm would stop, unable to reach that chain
at the top
> right corner. The algorithm has thus clearly not found
the proper length.
>
>
> Generally speaking, the algorithm's preference for
continuing in
> the same direction as possible is going to create
problems for you,
> as it will be unable to recognize where the molecule
changes
> direction even just due to a standard convex curve. And
as I showed
> above, counting pixels visited is not an accurate length
measure
> if any branches exist. Thirdly, you are not talking into
account
> that diagonal moves are further than vertical or
horizontal moves.
> --
> "Prevention is the daughter of intelligence."
> -- Sir
Walter Raleigh


Adrian,

This is very roughly coded but it gives 'an' answer. You
might find it helpful.

Ken

% Read in image and compress to 2D
ifs='c:\temp\mol.png';
im=imread(ifs);
im=sum(im,3);

% Show original image
figure(1);
imagesc(im);

% Blur and display
filt=fspecial('disk',2);
blurred_im=imfilter(im,filt,'replicate');

figure(2);
imagesc(blurred_im);
colorbar;

% Threshold the blurred version
threshold_im=zeros(size(im));
threshold_im(blurred_im>230)=1;

figure(3);
imagesc(threshold_im);

% And keep only the largest 'island' which is the molecule
molecule_im=bwareaopen(threshold_im,600);

figure(4);
clf;
imagesc(molecule_im);
hold on;

% Now get the perimeter
mol_properties=regionprops(bwlabel(molecule_im), ...
    'perimeter');
molecule_length=(mol_properties.Perimeter)/2

% And draw it to make sure we have the right thing
[p,l]=bwboundaries(molecule_im,4,'noholes');
boundary=p{1};
plot(boundary(:,2),boundary(:,1),'y');

Thank you so much for your responses and analysis of my
problem. I will try all of your suggestions and will post
anything I come up with.

Ken, thank you for the code - I just tried and it and seems
to be working very well. Your algorithm is very clean,
something I admire. I clearly need to learn more of the
advanced functions MATLAB has to offer.

Thanks again to everyone for your help,
Adrian Randall

"Ken Campbell" <campbeks@gmail.com> wrote in message
<g0d7jc$mtc$1@fred.mathworks.com>...
> roberson@ibd.nrc-cnrc.gc.ca (Walter Roberson) wrote in
> message <g0crs0$6al$1@canopus.cc.umanitoba.ca>...
> > In article <g0couj$np$1@fred.mathworks.com>, Adrian
> <ajr@med.unc.edu> wrote:
> >
> > >4) Scan the image from left to right and then "walk"
> along
> > >areas that are labeled as "1". If the program finds a
> pixel
> > >that is labeled "1", look toward the 8 surrounding
> pixels to
> > >find if there is another pixel labeled "1". If so,
> continue
> > >to look for pixels labeled "1" (e.g., keep "walking"
> along
> > >the molecule), and if possible, continue "walking" in the
> > >same direction as previously (to prevent the "walker"
> from
> > >going back to where it came from). Then, after walking
> along
> > >the molecule, label all of the points it walked on
> as "0" so
> > >the scanner won't walk there again.
> > >5) Output the length measurement as pixel and convert to
> nm
> > >by looking at the scale.
> >
> > Suppose you have
> >
> > 1011
> > 0100
> > 0010
> >
> > And you start at the top-left corner. You walk left to
> right so you
> > encounter the top-left 1. Presumably you increment your
> count to 1.
> > You examine the surrounding pixels and find
> > the 1 diagonally down and right. You try to continue left
> to right, but
> > there is no 1 there so you take the path that is open by
> going diagonally
> > down and right. However, as you were not able to continue
> in the same
> > direction as you had been going, you label the points you
> walked in
> > that stretch with a 0 so the scanner won't walk there
> again, leading to
> >
> > 0011
> > 0100
> > 0010
> >
> > You are at the second pixel, presumably you increment
> your count to 2.
> > Now you are headed diagonally down and right. You examine
> the
> > surrounding pixels and find the 1 diagonally up and right
> and
> > the 1 diagonally down and right. As per your algorithm,
> you prefer
> > to keep going in the same direction, so you take the
> diagonal
> > down and right. As you -were- able to continue in the
> same direction,
> > you are on a run and so do -not- mark the pixel you came
> from with a 0.
> >
> > You are at the third pixel, presumably you increment your
> count to 3.
> > You look around from that bottom pixel and find the only
> path open
> > to be the one up and left. However, as you were not able
> to continue
> > in the same direction as you had been going, you label
> the points
> > you walked in the run with a 0 won't walk there again,
> leading to
> >
> > 0011
> > 0000
> > 0000
> >
> > There is now an ambiguity in your algorithm: you had an
> available
> > path up-left a moment ago, but that pixel is now 0'd, and
> there
> > are now no pixels left in the neighbourhood that are non-
> 0. Do you
> > stop at this point with a pixel count of 3? My
> interpretation of
> > your algorithm is that you do -not- stop, because your
> instruction
> > was to do the search before the 0'ing.
> >
> > So you are now back at 2nd row 2nd column. Your algorithm
> didn't
> > say anything about not counting on backtracks, so
> presumably you
> > increment the count to 4. You look around and find the
> open path
> > to the up-right. That's a change of direction, so you 0
> out the pixels
> > you had in that run, which makes no change because they
> are already 0.
> >
> > You are now at the 1st row 3rd column, trying to head up
> and right.
> > Presumably you increment your counter to 5. You look
> around and
> > find the open path to the right. That's a change of
> direction
> > so you 0 out the pixels you had in that run, which makes
> the matrix
> >
> > 0001
> > 0000
> > 0000
> >
> > You are now at the 1st row 4th column, trying to head
> right.
> > Presumably you increment your counter to 6. You look
> around and find
> > no open paths. Presumably you stop there.
> >
> > Your recorded molecule length is now 6, but the longest
> connected
> > chain was 4 -- (1,1) -> (2,2) -> (1,3) -> (1,4), so
> arguably the
> > molecule length is 4 instead. Or, since there were two
> diagonal
> > moves and then 1 lateral move, 2*sqrt(2)+1, which is
> about 3.8.
> > Add 1 if you are counting from the edges instead of from
> the centers.
> >
> >
> > Now, suppose you had started with
> >
> > 101111
> > 010000
> > 001000
> > 000100
> >
> > then when you are examining (2,2) after coming from (1,1),
> > because your algorithm always continues in the same
> direction if it
> > can, it would continue on the downward diagonal to the
> edge, find
> > the upward-left step from the bottom as the only
> connection to the
> > bottom pixel, and would 0 out the entire diagonal,
> leaving you with
> >
> > 001111
> > 000000
> > 000000
> > 000000
> >
> > with you positioned at (3,3). There are no connections
> from there
> > so the algorithm would stop, unable to reach that chain
> at the top
> > right corner. The algorithm has thus clearly not found
> the proper length.
> >
> >
> > Generally speaking, the algorithm's preference for
> continuing in
> > the same direction as possible is going to create
> problems for you,
> > as it will be unable to recognize where the molecule
> changes
> > direction even just due to a standard convex curve. And
> as I showed
> > above, counting pixels visited is not an accurate length
> measure
> > if any branches exist. Thirdly, you are not talking into
> account
> > that diagonal moves are further than vertical or
> horizontal moves.
> > --
> > "Prevention is the daughter of intelligence."
> > -- Sir
> Walter Raleigh
>
>
> Adrian,
>
> This is very roughly coded but it gives 'an' answer. You
> might find it helpful.
>
> Ken
>
> % Read in image and compress to 2D
> ifs='c:\temp\mol.png';
> im=imread(ifs);
> im=sum(im,3);
>
> % Show original image
> figure(1);
> imagesc(im);
>
> % Blur and display
> filt=fspecial('disk',2);
> blurred_im=imfilter(im,filt,'replicate');
>
> figure(2);
> imagesc(blurred_im);
> colorbar;
>
> % Threshold the blurred version
> threshold_im=zeros(size(im));
> threshold_im(blurred_im>230)=1;
>
> figure(3);
> imagesc(threshold_im);
>
> % And keep only the largest 'island' which is the molecule
> molecule_im=bwareaopen(threshold_im,600);
>
> figure(4);
> clf;
> imagesc(molecule_im);
> hold on;
>
> % Now get the perimeter
> mol_properties=regionprops(bwlabel(molecule_im), ...
> 'perimeter');
> molecule_length=(mol_properties.Perimeter)/2
>
> % And draw it to make sure we have the right thing
> [p,l]=bwboundaries(molecule_im,4,'noholes');
> boundary=p{1};
> plot(boundary(:,2),boundary(:,1),'y');
>

Ken,

Is there anyway to modify your code to outline more than one
molecule? It works great for one molecule, but when I start
looking at more than one, the program breaks.

I'm planning on making a GUI for the program so the user can
change the thresholds as needed for an image. After trying
out your algorithm with various DNA images I have, the
shading is different between each one (because of the
inconsistent heavy-metal shadowing, I assume).

Thanks,
Adrian

<SNIP>

> > Adrian,
> >
> > This is very roughly coded but it gives 'an' answer.
You
> > might find it helpful.
> >
> > Ken
> >
> > % Read in image and compress to 2D
> > ifs='c:\temp\mol.png';
> > im=imread(ifs);
> > im=sum(im,3);
> >
> > % Show original image
> > figure(1);
> > imagesc(im);
> >
> > % Blur and display
> > filt=fspecial('disk',2);
> > blurred_im=imfilter(im,filt,'replicate');
> >
> > figure(2);
> > imagesc(blurred_im);
> > colorbar;
> >
> > % Threshold the blurred version
> > threshold_im=zeros(size(im));
> > threshold_im(blurred_im>230)=1;
> >
> > figure(3);
> > imagesc(threshold_im);
> >
> > % And keep only the largest 'island' which is the
molecule
> > molecule_im=bwareaopen(threshold_im,600);
> >
> > figure(4);
> > clf;
> > imagesc(molecule_im);
> > hold on;
> >
> > % Now get the perimeter
> > mol_properties=regionprops(bwlabel(molecule_im), ...
> > 'perimeter');
> > molecule_length=(mol_properties.Perimeter)/2
> >
> > % And draw it to make sure we have the right thing
> > [p,l]=bwboundaries(molecule_im,4,'noholes');
> > boundary=p{1};
> > plot(boundary(:,2),boundary(:,1),'y');
> >
>
> Ken,
>
> Is there anyway to modify your code to outline more than
one
> molecule? It works great for one molecule, but when I
start
> looking at more than one, the program breaks.
>
> I'm planning on making a GUI for the program so the user
can
> change the thresholds as needed for an image. After trying
> out your algorithm with various DNA images I have, the
> shading is different between each one (because of the
> inconsistent heavy-metal shadowing, I assume).
>
> Thanks,
> Adrian

You might be able to get round the inconsistent shadowing
by doing some background correction. Look at imtophat
and/or imbothat.

Getting the code to identify lots of molecules at the same
time is too complicated for the quick code I posted (your
original image only had one molecule after all), but the
general idea is the same. (1) Make the molecules stand out
by some blurring/background correction, (2) set the minimum
size of the molecule you want to analyze with the second
parameter to bwareaopen and (3) analyze the perimeters. You
might want to correct the 4/8 connected ness problem my
original code had.

I'm a bit too busy to do the whole thing for you :-) but it
shouldn't be too hard. Writing the GUI to do it
interactively will probably take you longer (it would
certainly take me longer) than getting the image-processing
part of the code to work once you understand the algorithm
and the underlying Matlab commands.

Ken