Help with identifiny change in curve shape automatically

Hi Everyone, I have a problem that should be pretty easy to solve. I am doing image analysis of a surface deformation and there are a few aspects of the surface that I want to analyze. The surface has a conical outer section and a parabolic bottom section. I want to extract the approximate points where the conical section meets the parabolic section. The following is a link to a picture of the curve (blue) and arrows indicating the information I am looking for http://www.flickr.com/photos/43987619@N04/6993746574/

I have tried working with the 1st and second derivative, but the results are not good. I wonder if there is a way I can post the actual curve points so you can try to play with them.

Thanks,

Casey

[EDIT: Used code markup so D is more easily pasted into MATLAB. - the cyclist]

D=[741  740  739  738  737  736  735  734  733  732  731  730  729  728  727  726  725  724  723  722  721  720  719  718  717  716  715  714  713  712  711  710  709  708  707  706  705  704  703  702  701  700  699  698  697  696  695  694  693  692  691  690  689  688  687  686  685  684  683  682  681  680  679  678  677  676  675  674  673  672  671  670  669  668  667  666  665  664  663  662  661  660  659  658  657  656  655  654  653  652  651  650  649  648  647  646  645  644  643  642  641  640  639  638  637  636  635  634  633  632  631  630  629  628  627  626  625  624  623  622  621  620  619  618  617  616  615  614  613  612  611  610  609  608  607  606  605  604  603  602  601  600  599  598  597  596  595  594  593  592  591  590  589  588  587  586  585  584  583  582  581  580  579  578  577  576  575  574  573  572  571  570  569  568  567  566  565  564  563  562  561  560  559  558  557  556  555  554  553  552  551  550  549  548  547  546  545  544  543  542  541  540  539  538  537  536  535  534  533  532  531  530  529  528  527  526  525  524  523  522  521  520  519  518  517  516  515  514  513  512  511  510  509  508  507  506  505  504  503  502  501  500  499  498  497  496  495  494  493  492  491  490  489  488  487  486  485  484  483  482  481  480  479  478  477  476  475  474  473  472  471  470  469  468  467  466  465  464  463  462  461  460  459  458  457  456  455  454  453  452  451  450  449  448  447  446  445  444  443  442  441  440  439  438  437  436  435  434  433  432  431  430  429  428  427  426  425  424  423  422  421  420  419  418  417  416  415  414  413  412  411  410  409  408  407  406  405  404  403  402  401  400  399  398  397  396  395  394  393  392  391  390  389  388  387  386  385  384  383  382  381  380  379  378  377  376  375  374  373  372  371  370  369  368  367  366  365  364  363  362  361  360  359  358  357  356  355  354  353  352  351  350  349  348  347  346  345  344  343  342  341  340  339  338  337  336  335  334  333  332  331  330  329  328  327  326  325  324  323  322  321  320  319  318  317  316  315  314  313  312  311  310  309  308  307  306  305  304  303  302  301  300  299  298  297  296  295  294  293  292  291  290  289  288  287  286  285  284  283  282  281  280  279  278  277  276  275  274  273  272  271  270  269  268  267  266  265  264  263  262  261  260  259  258  257  256  255  254  253  252  251  250  249  248  247  246  245  244  243  242  241  240  239  238  237  236  235  234  233  232  231  230  229  228  227  226  225  224  223  222  221  220  219  218  217  216  215  214  213  212  211  210  209  208  207  206  205  204  203  202  201  200  199  198  197  196  195  194  193  192  191  190  189  188  187  186  185  184  183  182  181  180  179  178  177  176  175  174  173  172  171  170  169  168  167  166  165  164  163  162  161  160  159  158  157  156  155  154  153  152  151  150  149  148  147  146  145  144  143  142  141  140  139  138  137  136  135  134  133  132  131  130  129  128  127  126  125  124  123  122  121  120  119  118  117  116  115  114  113  112  111  110  109  108  107  106  105  104  103  102  101  100  99  98  97  96  95  94  93  92  91  90  89  88  87  86  85  84  83  82  81  80  79  78  77  76  75  74  73  72  71  70  69  68  67  66  65  64  63  62  61  60  59  58  57  56  55  54  53  52  51  50  49  48  47  46  45  44  43  42  41  40  39  38  37  36  35  34  33  32  31  30  29  28  27  26  25  24  23  22  21  20; 78  77  78  79  79  79  79  80  80  80  80  80  80  81  81  81  81  81  81  81  82  82  82  82  83  83  83  83  83  83  83  83  83  84  84  84  84  84  84  84  84  84  85  85  85  85  85  85  85  85  85  85  85  85  85  85  85  85  85  85  85  85  85  85  86  86  86  86  86  86  86  86  86  86  86  86  86  86  86  86  86  86  86  86  86  86  86  86  86  86  86  86  86  86  86  86  86  86  86  86  86  86  86  86  86  86  86  86  86  86  86  86  86  86  86  86  86  86  86  86  86  86  86  86  86  86  86  86  86  86  86  86  86  86  86  86  86  86  86  86  86  86  86  86  86  86  86  86  86  86  86  86  86  86  86  86  86  86  86  86  87  87  87  87  87  87  87  87  87  88  88  88  88  89  89  90  90  90  91  91  92  92  93  94  94  95  95  96  96  97  97  98  98  99  99  100  100  101  101  102  102  103  103  104  104  105  105  106  106  107  107  107  108  108  109  109  109  110  110  111  111  112  112  113  113  114  114  115  115  115  116  116  117  117  118  118  119  119  119  120  120  121  121  122  122  123  123  124  124  125  125  125  126  126  127  127  128  128  128  129  129  130  130  130  131  131  132  132  133  133  134  134  135  135  136  136  137  137  137  138  138  139  139  139  140  140  140  141  141  141  142  142  142  143  143  143  144  144  144  145  145  145  146  146  146  147  147  148  148  149  149  150  151  152  153  155  156  158  160  162  164  166  168  169  171  172  174  176  177  178  180  181  182  183  184  185  187  188  189  190  191  192  193  194  195  195  196  197  197  198  198  199  199  200  200  200  200  200  201  201  201  201  201  201  201  201  201  201  201  201  201  201  201  201  201  201  201  200  200  200  200  199  199  198  198  197  196  196  195  194  193  192  191  191  190  189  188  187  186  185  184  183  182  180  179  178  176  175  174  172  171  169  168  166  164  163  161  159  158  156  155  153  152  151  151  150  150  149  149  148  148  148  147  147  146  146  146  146  145  145  144  144  144  143  143  143  143  142  142  142  141  141  140  140  140  139  139  139  138  138  137  137  136  136  136  135  135  134  134  134  133  133  132  132  132  131  131  130  130  129  129  128  128  127  127  127  126  126  125  125  124  124  123  123  122  122  122  121  121  120  120  119  119  119  118  118  118  117  117  117  117  116  115  115  114  114  113  113  112  112  111  111  110  110  109  109  109  108  108  107  107  106  106  105  105  104  104  104  103  102  102  102  101  101  100  100  100  99  99  98  98  97  97  96  96  95  95  94  93  93  92  92  91  91  91  90  90  90  89  89  89  88  88  88  87  87  87  87  87  87  87  87  87  87  86  86  86  86  86  86  87  87  87  87  87  86  86  86  86  86  86  86  86  86  86  86  86  86  86  87  87  87  87  87  87  87  87  87  87  87  87  87  87  87  87  87  87  87  87  87  87  87  87  87  87  87  87  87  87  87  87  87  87  87  87  88  88  88  88  88  88  88  88  88  88  88  88  88  88  88  88  88  88  88  88  88  88  88  88  88  88  89  89  89  89  89  89  89  89  89  89  89  89  89  89  89  89  89  89  89  89  89  89  89  89  89  89  89  89  89  89  89  89  89  89  89  89  89  89  89  89  89  89  89  89  89  89  89  88  88  88  88];

[Image Analyst]

Is this ideal data, or will there be some amount of noise on that signal that will make it more difficult to identify the "kink" in it?

[Casey]

This is pretty typical data. Any data that is worse than this I will likely have to do by hand. There are some times when the conical part of the surface slides down and there is a small bump separating the conical portion from the parabolic portion. I can try and find a set of data to show this, but it may take a while.

Answer 1

[Richard Brown]

What I think you're looking for are effectively the points of high curvature. If you have access to the curve fitting toolbox then you can attack it with an aggressive smoothing spline to make the derivatives smooth, and then compute the curvature. The points of interest are then the local maxima of the curvature.

Here's some code that hopefully does the trick. I'll start with the matrix D from above

x = fliplr(D(1,:));
dx = 1;
ynoisy = fliplr(D(2,:));

Now we'll smooth it - you may need to tinker with the smoothing parameter - closer to zero is smoother.

p = 1e-4;
fn = csaps(x, ynoisy, p);
y = ppval(fn,x);

Take a quick look to make sure it's ok

plot(x, ynoisy, 'rx', x, y, 'b')

Now we'll compute first and second derivatives using centred finite differences

yp = nan(size(y));
ypp = nan(size(y));
yp(2:end-1) = (y(3:end) - y(1:end-2))/(2*dx);
ypp(2:end-1) = (y(3:end) + y(1:end-2) - 2*y(2:end-1)) / (dx^2);

Curvature is then a simple computation

k = abs(ypp) ./ (1 + yp.^2).^1.5

If you plot the curves together, you can see that the peaks of k correspond to the turning points that you are interested in

plotyy(x, ynoisy, x, k)

If you hvae the signal processing toolbox too, then you'll have access to findpeaks too. If not, there's probably something on the file exchange. At any rate, what you're interested in is the five highest peaks

[pks, iPeaks] = findpeaks(k);
[~, iSorted] = sort(pks, 'descend');
iPeaks = iPeaks(iSorted(1:5));

Now plot them to make sure we've found them correctly

plot(x, ynoisy, 'b', x(iPeaks), ynoisy(iPeaks), 'ro')

[Casey]

This works great! Thanks!! Can you explain why you found the derivatives that way? I think that may be why my attempts were failing before.

[Casey]

Now, if the data is noisier, I would just increase p, right?

[Richard Brown]

The derivatives are just standard central finite difference approximations (you should be able to find them on Wikipedia). I prefer to use a central difference for the first derivative as with a small dx it is more accurate.

The critical thing is that you get your data smooth enough, otherwise any derivatives you compute will be garbage. The smoothing parameter for csaps is pretty sensitive, and *decreasing* it makes the curve smoother. If you're having trouble tuning it you could try using spaps instead, which has pretty similar functionality. You might also notice that the peak in the middle of your plot turns into two peaks sometimes, depending on your choice of parameter - so you might need to tweak the last bit too!

[Casey]

Thanks Richard, this is working very well for the most part. Some of my videos have more noise and it makes it hard to find the correct changes in curvature. I was wondering if there is a way to fit the curve to two straight lines separated by a parabola. Then I can define the points I am looking for as the points that separate the straight lines from the parabola. If you would like to PM back and forth let me know and I will try and send you my email. I really appreciate all of your help!

[Richard Brown]

Sure, feel free to drop me an email