How to redistribute pixels?
Show older comments
Posted is a picture of a cross section of a vein. I have transformed the red interior line into a circle and need to have the pixels on the black line surround the new circular red line, both red and black lines having the same centers. The transformed black line should not be a perfect circle. As you can see in the picture, the distance between the black and red lines differ and once the black line is transformed, it should maintain the same varying distance between it and the red line as before. My question is how can I do this?
Answers (2)
Image Analyst
on 27 Jan 2013
1 vote
It's not too hard if you think about it systematically. Just take it one step at a time. The key is to convert the lines into masks and use bwdist().
- use poly2mask on the black line to get a mask
- use poly2mask on the red line to get a smaller mask
- subtract the two masks to get the "in between" region.
- call bwdist() to get the Euclidean Distance Transform (EDT). This is the distance from the region to the nearest edge, so twice that is the thickness of the region at that point. So multiply the EDT by 2.
- call bwmorph() to get the skeleton. This is the centerline of the region.
- Multiply the skeleton by the EDT. This will give you the thickness of the region.
- bwdist will also give you the pixel on the red line that is closest to the black line.
- Now you know the radius and center of the red circle. Radius is just perimeter/(2*pi). So you just run around the black circle, and get the red pixel that is closest to it.
- You know the position of each red pixel - it's index on the circle. So you plot a red point there, and then plot a black point at that red radius plus 2*EDT for that black pixel.
4 Comments
Tony D
on 28 Jan 2013
Image Analyst
on 29 Jan 2013
(4) The region thickness is twice the EDT value at the center. Let's say you had a region 20 pixels wide. Well, at the center it will be 10 pixels to either side because if you're at the center line it has to go only half the thickness to get outside the region. So you need to double the value at the center to get the full thickness.
6) You multiply by the skeleton to give you the value along the centerline of the region. The EDT you get in (4) is a thick image. It's like a wide swatch the same size as the region in between your two curves. And the values of the EDT go from 0 or 1 at the edge (where a pixel inside the swath does not have far to go to get outside) all the way to half the thickness at the centerline. If you're 3 pixels away from the background, the EDT value is 3, if you're 4 pixels away it's 5, if you're 5 pixels away it's 5, and so on until you're the radius (half thickness) away from the edge. If you proceed toward the other side, it starts going down again until it again reaches 0 at the other edge (boundary).
Think of it like this. You're taking an inflatable ball and running it along the region between the two lines. The ball is always inflated as much as it can go and still be a circle. Sometimes, when the thickness is large, you can fit a large ball, and sometimes, when the thickness is small, you'll have a small ball. The EDT along the centerline (skeleton, spine) is the radius of the ball at that point.
Does that explain it better?
Tony D
on 30 Jan 2013
Tony D
on 30 Jan 2013
Matt J
on 4 Feb 2013
So now how can I use MATLAB to go through all of the pixels on a boundary of the matrix and find the distance to the nearest pixel on the opposing line?
If the rows of A are the pixel coordinates on the black line and B the red line, then you can use
as follows
s=ipdm(A,B,'Subset', 'NearestNeighbor','Result','Structure');
Then s.columnindex will be the indices in B of the closest points on the red line. The corresponding distances will be in s.distance.
When run with this specific syntax, the column index is what will be most useful. It tells you that
B(short_dist.columnindex(i), :)
is the closest point on the red line to A(i,:). And if you want the vector between the two points, just do
B(short_dist.columnindex(i), :) - A(i,:)
for all i=1...1856. Or, in vectorized form, simply
B(short_dist.columnindex, :) - A
42 Comments
I implemented your idea into my code and it almost works. It has some difficulty calculating the new red line around the right side of the wall where it is folded up. Here is how I am calculating the transform so far:
- use the ipdm function to calculate the vectors between the red and black lines pre-transformation using your previous comment
- calculate the perimeter of the black line [units=pixels]
- convert that perimeter into a radius
- multiply that radius by the pixel resolution to convert units into microns
- using that radius, draw a circle that is centered at the origin (has the same number of points that the black line had in pixels)
- add the vectors from the ipdm function to the circle to get the red line.
As you will see in the picture, the boundaries are not correct. Any ideas on what when wrong?
I' a bit confused. In your originally posted question, you say the red line is transformed into a perfect circle and the black line transforms into a non-circle derived from the red line. Here, however, you seem to have done the opposite.
I suspect you will get what you're looking for if you just interchange the roles of A and B in all steps, including the use of IPDM, i.e., you should be finding the distances/vectors of red-to-black instead of black-to-red.
Matt J
on 5 Feb 2013
Tony DeStephens commented:
Yes, originally I wanted to transform the red line into a perfect circle, now I want to instead transform the black line into a perfect circle and the red line which is a non-circle gets derived from the black line. I apologize for the confusion.
So if I now wanted to instead keep the black line a circle I would just re-arrange the roles of A and B? Something like this:
s=ipdm(B,A,'Subset', 'NearestNeighbor','Result','Structure');
vector = A(short_dist.columnindex(i), :) - B(i,:);
A_new = equations for cirlce;
B_new = A_new + vector;
Is this correct?
So if I now wanted to instead keep the black line a circle I would just re-arrange the roles of A and B?
No. I only suggested interchanging A and B because it sounded like you had the red line and black line mixed up. If you're now interested in transforming the black line into the circle, that doesn't apply.
6.add the vectors from the ipdm function to the circle to get the red line.
I think this step is the real problem. It sounds like you want to be stepping towards the center of the black circle by the length of the IPDM vectors. In particular, let v be the vector given by ipdm, let A_new be the corresponding point on the black circle, and let C be the position of the center of the black circle. Then I think you want,
B_new = A_new + norm(v)*(C-A_new)/norm(C-A_new)
Matt J
on 6 Feb 2013
I had v = B(short_dist.columnindex, :) - A
No, you need to be looping over all A_new individually.
Incidentally, you should be replying to me in comments so that people know whose answer you're addressing.
Tony D
on 11 Feb 2013
Your code worked great, thanks for the help!
To check that the transform was done properly, I am calculating the area between the two lines both before and after the transformation and comparing them. I am using the poly2mask function to create a binary image from the post transform coordinates, but only 1/4 of the image is showing up. I am using the following commands where BW is the binary image of the ROI pre-transformation:
BW_A_new = poly2mask(A_new_xunit, A_new_yunit, length(BW(:,1)), length(BW(1,:)));
BW_B_new = poly2mask(B_new_xunit, B_new_yunit, length(BW(:,1)), length(BW(1,:)));
in_between_region_post = BW_A_new - BW_B_new;
figure (4), imshow(in_between_region_post)
Am I calling poly2mask incorrectly? How do I get the whole image to so up?
Also, I am using the imfreehand function to draw the lines on the picture, but I am going to have to draw a 3rd line that outlines the boundary of the inside diameter (which is just inside the red line). In order to do this I will have to be constantly stopping, zooming in and moving the image. So the question is how can I stop and go while outlining an image using imfreehand?
Matt J
on 11 Feb 2013
I think it's because you're using x,y coordinates with an origin somewhere in the middle of the image. They need to be in matrix coordinates, if I'm not mistaken.
I'm a bit skeptical though that the area between the curves will be preserved (I could be wrong, though). Was that the goal all along? You should have said so in your original post.
Tony D
on 12 Feb 2013
how do I make a transformation from Cartesian coords to matrix coords?
You first convert Cartesian coordinates to pixel length units (instead of millimeters or whatever). Then you translate to a new origin. The matrix coordinate origin is in the top left corner.
If the distances between the pixels are the same before and after the transform then my thought was that the area change would not change.
Don't think so. Consider the simple case of 2 concentric circles with radii a>b. The area between them is
InitialArea = pi*(a^2-b^2)
Now transform the circles simply by increasing both radii by the same amount, x, so that the distance between them is preserved at a-b. The new area is
NewArea = pi*((a+x)^2-(b+x)^2)
= InitialArea+2*pi*(a-b)*x
Clearly the distance between the curves is invariant to x but the area between them is not.
Tony D
on 12 Feb 2013
Matt J
on 13 Feb 2013
When I did this and used imshow to bring up the image nothing showed up. The figure was completely black and the ROI was gone.
We can't see your calculations, so we can't see the relationship between your current coordinates and the one that poly2mask works with. By experimenting with with poly2mask and reading the examples in its documentation, you should be able to understand the coordinates it uses.
I see what you are saying, but what if the perimeter of the black line was preserved, so that it was the same for the pre and post transform as well as preserving the distances between the red and black line?
It might be true if the closest point from points on the A curve to the B curve is a 1-1 mapping. Not sure how to confirm that with rigorous calculus, though. It doesn't seem to be true when the mapping is not 1-1 however.
Suppose for example your red curve was 1 side of a unit square and your black curve were the remaining 3 sides of the square. In this case, the tips of the red edge are the nearest points to every point on the black edge intersecting them, so the nearest-point mapping is not 1-1. Now think about how the initial area (equal to 1) changes if you transform the black curve by straightening it into a line segment with the same initial total length = 3. If you construct a new red border based on distance along rays perpendicular from this new black line, you will end up with a trapezoid of area 2.
Also, how can I stop and go while outlining an image using imfreehand?
Don't think you can do that, but you can call imfreehand with 'Closed'=false. Then it will just draw an open curve. You can then zoom and call imfreehand again and draw a new segment of the curve, etc... With just a little bit of coding, you can accumulate all the points on the different curve segments.
Tony D
on 19 Feb 2013
Can you explain how you got this equation that calculates the new red line?
The length of the new difference vector (B_new - A_new) will be norm(v)=norm(B_old-A_old), preserving the previous distance. The direction of the new difference vector will be (C-A_new)/norm(C-A_new) which points radially inward from the boundary of the new black circle.
I have now implemented this to my code, but is there a way that I can have a newly drawn open line to automatically connect the previously drawn open line?
Not that I know of, but you could consider using impoly instead of imfreehand. With impoly, you can freely zoom in, add new vertices, and bend the boundary of the polygon as finely as you like.
Then maybe add that vector to the transformed point on the red line?
Thought you tried that before and it wasn't what you wanted.
Tony D
on 26 Feb 2013
Would this be a better method in preserving the distance between the two lines than finding the nearest neighbor
No one but you knows what "better" means (because it depends on what you're trying to achieve with all this).
one of the goals for this transfiguration is to have a 1-1 mapping.
It would really be better if you described all of the goals of the transformation. As you can see, leaving design requirements unspoken leads to proposals that don't satisfy them.
So I was thinking that instead of going pixel by pixel on the black line and finding the nearest pixel neighbor on the red line, instead find the red line pixel that is orthogonal to the black line pixel and tracking that distance instead.
That wouldn't guarantee a 1-1 mapping, I'm afraid. If either curve bends sharply enough then a perpendicular emanating from certain points on the outer black curve may not intersect the interior red curve at all (imagine concentric squares). Conversely, a perpendicular radiating from the interior curve to the outer curve would skip over points at the corners.
I think you need to be analyzing this more generally. Suppose we parametrize the given red and black curves, (e.g., using splines) as
r_black(t) = [x_black(t), y_black(t)]
r_red(t) = [x_red(t), y_red(t)]
Any such parametrization defines a 1-1 mapping between r_red and r_black as long as each curve is swept out completely as t increases monotonically from 0 to 1.
Similarly write down a parametric formula for the target black circle
r_circ(t) = [x_circ(t), y_circ(t)]
Then there is an infinite set of distance-preserving new red cuves r(t) given by
r(t) = r_circ(t)+ norm(r_red(t) - r_black(t))* v(t)
where v(t) is any vector satisfying norm(v(t))=1 at all t. The set is infinite because there are infinite possibilities for v(t) and also because there are an infinite variety of ways to parametrize r_red, r_black, and r_circ as functions of t.
The thing to do therefore is to specify what properties that you want r(t) to have in addition to distance preservation (all choices satisfy that) so that the choices for some of these things can be narrowed down.
Tony D
on 5 Mar 2013
So if the red and black lines were drawn as monotonic splines, then 1-1 mapping will be achieved?
I don't know what you mean by a "monotonic spline". If you have 2 space curves r1(t) and r2(t) which are each 1-1 functions of t, then likewise there is a 1-1 relationship between the points on the curves, given by
r2(t(r1))
In other words, a given coordinate r1 maps to some unique t which maps to a unique r2 and vice versa.
Other than that I want the region between the lines pre-transformation to be the same as post-transformation.
What does it mean for 2 regions to be "the same"? Obviously the regions will have a different shape post-transformation. You mean their areas should be the same?
Tony D
on 5 Mar 2013
Matt J
on 6 Mar 2013
So right now the red and black lines are hand drawn, but your saying that they need to be defined by parametric equations each being a function of t?
Yes, but for example, if you used impoly instead of imfreehand, you know that the curve is piecewise linear. It would be very easy to write down a piecewise linear r(t) once the vertices are obtained. Similarly, you could use the SPLINE command to get a cubic spline fit of the vertices
Yes, I want the areas to be the same. sorry for the confusion.
This is starting to sound doubtfully difficult. In order for the transformation to be both length and area preserving, it would have to probably consist of piece-wise rigid rotation/translations that would map small sections of the initial black curve to the final black circle. If you applied such a transformation, I can see the red curve getting transformed in ways you don't like, e.g., bending in ways that make it self-intersecting...
What exactly is the practical purpose of all this again? If we succeed at finding this transformation, what does it then allow you to do?
Matt J
on 9 Mar 2013
I basically want to track the radial position of the cells in between those two lines after the transformation.
So basically, the whole idea all along was to transform the 2 curves into polar curves r(theta) where r is a 1-1 function of theta? Does the transformed black curve even have to be a circle?
Tony D
on 9 Mar 2013
Does the transformed black curve even have to be a circle? Yes
I was looking for a little more elaboration here. It doesn't sound like it needs to be a circle, if all you want to be able to do (post-transformation) is to draw a radial line through each cell and see where that line intersects the red and black curves.
Aside from that, I'd like to step back for a moment and try to reverse engineer things. Suppose we start with the transformed red and black curves that you're looking for. Suppose I draw a radial line to the black circle and observe the chord running between its intersection with the red and black curves. Now, suppose that I invert the transformation, taking us back to the original red and black curves, but I also apply that inverse transformation to the chord (giving me a "pre-chord") and observe how it changes shape.
My question(s): what hard requirements do you have on that pre-chord? You've said you want the pre-chord and the chord to have the same length. However, does the pre-chord still have to be a straight line? Should the pre-chord intersect the pre-red and pre-black curves at only its 2 end points or is it allowed to cross them at other places as well?
Obviously, the more conditions you have on how chords/pre-chords map to one another, the harder it will be to find the desired transformation (if one even exists), so you should set only the requirements that you really and truly need. You should also explain why you think you need those conditions, so I can challenge that reasoning (some of the things you want may turn out to be impossible) and maybe find more feasible alternatives.
Tony D
on 11 Mar 2013
When the vein has blood is flowing through it, the black line, which represents a boundary in the vein wall, is circular
Do you know for a fact that the "distance" between the red and black curves remains the same in physical reality? If you know that for a fact, you ought to have some basis, based on biology, of defining that distance as well as the point on the red curve that a given point on the black curve is partnered with. We've been lacking those definitions, so far. Obviously, the nearest-point mapping from the black curve to the red curve does not define the partnering. We tried that...
To my inexpert eye, I can see that the total area between the 2 curves would remain the same assuming the material sandwiched in between the 2 curves (cells?) has a fixed density, but otherwise no other biomechanical laws about how the red boundary should transform are apparent to me. We need to know those laws.
Yes, but only for later transformations. Right now, for the first of many transforms, these laws are not taken into account.
I think you might be missing my point. I'm suggesting that not taking at least some of those laws into account is making things harder than including them. Even if you just want to preserve area and "distance" (whatever that is) between the two boundaries, the reason those things are preserved in real life is the complex result of fluid mechanical phenomena like applied pressure and non-compressibility. Writing down the equations for those gives one a better chance of deriving the proper transformation.
I think we're miles away from accomplishing this with an ad hoc geometric transformation. Rigidly rotating the black boundary piece by piece to align with the desired circle doesn't seem to me like it's going to work. If you were to do that in real life, it could stretch (or tear) the red boundary segment opposite it and that would change area.
Tony D
on 18 Mar 2013
I have a feeling you still imagine the region between the curves as consisting of trapezoidal sections. For trapezoidal partitions to be possible, the black curve would have to be approximately parallel to the nearby red curve everywhere. But I see parts of the curves in your picture where this is distinctly not the case.
When referring to stretching, I was talking about the earlier idea of rotating different sections of the curves, though in my mind, the sections were quadrilaterals, not specifically trapezoids. To bend the black curve into a circle, you would have to rotate a quadrilateral section about the point on the black curve where neighboring sections meet, but that would exert stress on the point on the red curve where the same two sections meet.
I'm afraid I've lost faith in IPDM as the way to define corresponding points on the red and black curves. As I said earlier, I'm not confident that it maps point A on the black curve 1-1 to a point B on the red curve. I'm also not sure that all such chords AB do not cross. If chords do cross, those chords cannot map to distinct chords after the black curve is bent into a circle, contrary to your goals. You're hoping, I think, that all the chords map to distinct radial lines under that transformation...
All I need is the unpressurized circular geometry,
The "unpressurized circular geometry" is new terminology. How does the black boundary become circular if not by the application of some sort of pressure? And I thought you said earlier that it was preciesly the bloodflow that makes it circular in the first place, so what do you mean by "how the geometry changes in the blood flow environment?
Can you see any way in which I can do the transform that I want using the given specifications?
Well, here's another idea. You're already using imfreehand or similar to trace the boundary of the curves. It seems to me that it would be essentially the same manual effort to use IMLINE or similar to draw lots of closely spaced line segments connecting the black curve to the red curve (essentially the pre-chords that we've been talking about). You would do this according to your own visual idea of which points on the black and red curves correspond, as opposed to using IPDM.
Once you've done that, you map all the points on the black curve to the desired circle. You would probably do this so that the distances between neighboring black points are all the same post-transformation as they were pre-transformation. You then use the lengths of the pre-chords to construct post-chords of the same length, except that the post chords are all directed inward toward the center of the black circle. The inner endpoints of the post-chords define the points on the transformed red curve and you're done.
Assuming you pack the pre-chords densely enough, you should get a pretty good approx of the curve you're looking for. Assuming, further, that a transform preserving both chord length and area does exist, I would expect this manual approach to do a fair bit better than the 8% area agreement you got via IPDM. If you don't get good area preservation, it would suggest to me that the requirements of chord length preservation and area preservation are in conflict. In either case, I wonder if you should really be caring about preserving distance and area so precisely. It sounds to me like you're planning to use this anyway as nothing more than an ad hoc initial approximation to something that's going to be refined later using something more principled (FEA?).
Tony D
on 25 Mar 2013
Matt J
on 25 Mar 2013
I don't know what it means to "multiply a curve" by something. If you mean to scale the chord lengths by that, then perhaps, but then chord-length won't be preserved.
I ran the code with a higher resolution photo (10710 x 6757) of a different case. Below are the binary images of the vein before and after the remapping. As you can see, the distances between the curves are not even close to matching. Any ideas why?
Also, can you think of a method that tests the output of the code to make sure that it is doing what it is suppose to be doing? The method that I have been using so far is to just take the original geometry and pick two points that seem like they are the furthest distance apart and find the magnitude of the distance between them. Then I see if there are a similar set of points of the remapped geometry
It's not clear how you did the remapping that produced these 2 most recent images. Did you use the approach I last suggested on Mar. 19?
In any case, you haven't displayed the curves with aspect ratio = 1 (e.g. with the command 'axis equal' or 'axis image'), so it's hard to judge how distances changed under the remapping. With Aspect Ratio = 1, the circle would look like a circle and not an ellipse.
I used the process that you recommended on Feb 4th, which used IPDM. I plotted the curves and then used 'axis equal' on both of the plots. I still am not sure if the IPDM worked the way that I wanted it to, but I do not know for sure and that is just by "eying it". Can you think of a way to track a pair of points pre transform and see how they mapped post transform?
Can you think of a way to track a pair of points pre transform and see how they mapped post transform?
Well, if you perform the transform on a single pair of points at the command line, you should be able to find them on the before/after plots just by looking at the coordinates that the transformation gives you.
For now, compressibility effects are assumed to be non existent. All I want for this first transformation is to preserve the distance and area between the two curves, make the black line circular,
I was just looking back at March 11. It seems like you were implying here that incompressibility implies the preservation of curve separation "distance", however that's supposed to be defined (I still think you haven't supplied a definition). Incompressibility obviously implies that the total area does not change (as long as the density of the material between the curves is uniform), but there is no clear distance preservation rule that this implies, as far as I can see, even combined with the fact that the outer curve's perimeter does not change.
As a counter-example, consider the case of 2 concentric squares with corresponding edges parallel. Let the outer square have area=2 and the inner square have area=1. Now transform by rotating the inner square 45 degrees so that its corners are now touching the centers of the sides of the outer square. It is simple to verify that the total area separating the squares has not changed, and obviously the outer square's perimeter hasn't changed. The transformation is also 1-1.
However, separation distances between the curves have obviously changed. For one thing, the minimum distance between the curves as now zero, i.e., where the inner square's corners touch the outer square. Corresponding points before and after the transformation, have also changed in separation distance, regardless of how you define that correspondence. For example, let's suppose that point A on the outer square corresponds to point B on the inner square (pre-transformation) if they have the same polar coordinate theta. Thus, initially, the corners of the outer square correspond with the corners of the inner square. The distance from the outer to inner corners changes from 0.2929 to 0.7071 after the rotation, however.
Conclusion: Non-compressibility plus constant outer perimeter doesn't imply separation distance preservation. So why are you pursuing it?
Tony D
on 10 Jun 2013
I was looking for a more efficient way than visual inspection. Maybe plotting individual points? I would just use the IPDM output to track the points?
It's not really clear to me what you're trying to "track". Are you just trying to verify that a pair of points A and B have the same separation distance after transformation? If so, you could just transform an individual pair of points and use the norm() function to verify that distance has remained the same.
This transformation would be similar to taking a wet spaghetti noodle and moving it (without stretching) to make a circle.
Okay. But what I was saying at that time is that the ability of wet spaghetti noodle to behave so malleably, and yet not dissipate completely like a gas, is probably still due to mechanical laws that aren't trivial. Otherwise, other things (e.g., dry spaghetti noodles and gas) would be able to do so as well. If we knew what those mechanical laws were, we might have better luck at identifying the "wet spaghetti noodle transform" in mathematical form :-)
It also seems to me that wet spaghetti must have some ability to stretch/compress to allow it to reshape as it can, and that lengths/distances aren't perfectly preserved when it does so. For example, if you lay a spaghetti noodle out straight on a table, the length of the noodle is the same, no matter along which edge of the noodle you measure the length. However, if you then join the tips and reshape the noodle into a circle, then obviously the outer edge has a longer circumference than the inner edge. So, obviously lengths/distances have changed somewhere.
Categories
Find more on Get Started with Computer Vision Toolbox in Help Center and File Exchange
Products
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!
