[R] Idea/package to "linearize a curve" along the diagonal?
Emmanuel Levy
emmanuel.levy at gmail.com
Tue Mar 13 19:29:38 CET 2012
Dear David and Jeff,
> Only if you were going apply some sort of transformation that did not extend globally
Exactly, this is why the LPCM package is great, as it assigns points
to parts of a curve.
I think I pretty much got what I need - it is not perfect yet but it
should be enough to give you an idea of what I was trying to achieve.
All the best,
Emmanuel
### Example ###
library(LPCM)
tmp=rnorm(2000)
X.1 = 5+tmp
Y.1 = 5+ (5*tmp+rnorm(2000))
tmp=rnorm(1000)
X.2 = 9+tmp
Y.2 = 40+ (1.5*tmp+rnorm(1000))
X.3 = 7+ 0.5*runif(500)
Y.3 = 15+20*runif(500)
Y = c(X.1,X.2,X.3)
X = c(Y.1,Y.2,Y.3)
lpc1 = lpc(cbind(X,Y), scaled=FALSE, h=c(1,1) , control=lpc.control(
boundary=1))
my.proj = lpc.spline(lpc1, project=TRUE, optimize=TRUE)
data = cbind( dist= my.proj$closest.pi, X1=lpc1$data[,1],
Y1=lpc1$data[,2], Xo=my.proj$closest.coords[,1],
Yo=my.proj$closest.coord[,2])
transfoData = matrix(apply(data, 1, function(x) { return( transfo(
(5+x[1])/10,x[2],x[3],x[4],x[5]))}), ncol=2, byrow=TRUE)
plot(transfoData) ## This shows the result I'm looking for, not
perfect yet but it gives an idea.
###
### Moves a point from it's position to the new "normalized" position
###
transfo = function(dist, X1, Y1, X0, Y0) {
# First, the point needs to be rotated
trans=newCoord(X1, Y1, X0, Y0) ;
Xnew=X1+trans[1]
Ynew=Y1+trans[2]
# second it is taken on the diagonal.
Xfinal=dist
Yfinal=dist
X.TransToDiag = Xfinal-X0
Y.TransToDiag = Yfinal-Y0
return( c(Xnew+X.TransToDiag, Ynew+Y.TransToDiag))
}
## Rotates a point X1,Y1 relative to Xo,Yo
## The new point is either at 3pi/4 or 7pi/4 i.e., 90 degrees left or
## right of the diagonal.
##
newCoord = function(X1,Y1, Xo=0, Yo=0){
# First calculates the coordinates of the point relative to Xo,Yo
Xr = X1-Xo
Yr = Y1-Yo
# Now calculates the new coordinates,
# i.e.,
# if V is the vector defined from Xo,Yo to X1,Y1,
# the new coordinates are such that Xf, Yf are at angle TETA
# by default TETA=3*pi/4 or 135 degrees
To = atan2(Yr,Xr)
# XXX This is not perfect but will do the job for now
if(Yr > Xr){
TETA=3*pi/4
} else {
TETA=7*pi/4
}
Xn = Xr * (cos(TETA)/cos(To))
Yn = Yr * (sin(TETA)/sin(To))
# Xn, Yn are the new coordinates relative to Xo, Yo
# However for the translation I need absolute coordinates
# These are given by Xo + Xn and Y0 + Yn
Xabs = Xo+Xn
Yabs = Yo+Yn
## the translation that need to be applied to X1 and Y1 are thus:
Xtrans = Xabs-X1
Ytrans = Yabs-Y1
return(c(Xtrans,Ytrans))
}
On 12 March 2012 20:58, David Winsemius <dwinsemius at comcast.net> wrote:
>
> On Mar 12, 2012, at 3:07 PM, Emmanuel Levy wrote:
>
>> Hi Jeff,
>>
>> Thanks for your reply and the example.
>>
>> I'm not sure if it could be applied to the problem I'm facing though,
>> for two reasons:
>>
>> (i) my understanding is that the inverse will associate a new Y
>> coordinate given an absolute X coordinate. However, in the case I'm
>> working on, the transformation that has to be applied depends on X
>> *and* on its position relative to the *normal* of the fitted curve.
>> This means, for instance, that both X and Y will change after
>> transformation.
>>
>> (ii) the fitted curve can be described by a spline, but I'm not sure
>> if inverse of such models can be inferred automatically (I don't know
>> anything about that).
>>
>> The procedure I envision is the following: treat the curve "segment by
>> segment", apply rotation+translation to each segment to bring it on
>> the
>> diagonal,
>
>
> That makes sense. Although the way I am imagining it would be to do a
> condition (on x) shift.
>
>
>> and apply the same transformation to all points
>> corresponding to the same segment (i.e., these are the points that are
>> close and within the "normal" area covered by the segment).
>>
>> Does this make sense?
>
>
> The first part sort of makes sense to me... maybe. You are thinking of some
> sort of local transformation that converts a curve to a straight line by
> rotation or deformation. Seems like a problem of finding a transformation of
> a scalar field. But you then want it extended outward to affect the regions
> at some distance from the curve. That's where might break down or at least
> becomes non-trivial. Because a region at a distance could be "in the sights"
> of the "normal" vector to the curve (not in the statistical sense but in the
> vector-field sense) of more than one segment of the curve. Only if you were
> going apply some sort of transformation that did not extend globally would
> you be able to do anything other than a y|x (y conditional on x) shift or
> expansion contraction
>
>> All the best,
>>
>> Emmanuel
>>
>>
>> On 12 March 2012 02:15, Jeff Newmiller <jdnewmil at dcn.davis.ca.us> wrote:
>>>
>>> It is possible that I do not see what you mean, but it seems like the
>>> following code does what you want. The general version of this is likely to
>>> be rather more difficult to do, but in theory the inverse function seems
>>> like what you are trying to accomplish.
>>>
>>> x <- 1:20
>>> y <- x^2 + rnorm(20)
>>>
>>> y.lm <- lm( y ~ I(x^2) + x )
>>> plot( x, y )
>>> lines( x, predict( y.lm ) )
>>>
>>> qa <- coef(y.lm)["I(x^2)"]
>>> qb <- coef(y.lm)["x"]
>>> qc <- coef(y.lm)["(Intercept)"]
>>>
>>> # define inverse of known model
>>> f1 <- function( y ) { ( sqrt( 4*qa*( y -qc ) + qb^2 ) - qb ) / ( 2*qa ) }
>>> f2 <- function( y ) { -( sqrt( 4*qa*( y -qc ) + qb^2 ) + qb ) / ( 2*qa )
>>> }
>>>
>>> plot( x, f1( y ) )
>>>
>>>
>>>
>>> ---------------------------------------------------------------------------
>>> Jeff Newmiller The ..... ..... Go
>>> Live...
>>> DCN:<jdnewmil at dcn.davis.ca.us> Basics: ##.#. ##.#. Live
>>> Go...
>>> Live: OO#.. Dead: OO#.. Playing
>>> Research Engineer (Solar/Batteries O.O#. #.O#. with
>>> /Software/Embedded Controllers) .OO#. .OO#.
>>> rocks...1k
>>>
>>> ---------------------------------------------------------------------------
>>> Sent from my phone. Please excuse my brevity.
>>>
>>>
>>>
>>> Emmanuel Levy <emmanuel.levy at gmail.com> wrote:
>>>
>>>> Dear Jeff,
>>>>
>>>> I'm sorry but I do not see what you mean. It'd be great if you could
>>>> let me know in more details what you mean whenever you can.
>>>>
>>>> Thanks,
>>>>
>>>> Emmanuel
>>>>
>>>>
>>>> On 12 March 2012 00:07, Jeff Newmiller <jdnewmil at dcn.davis.ca.us>
>>>> wrote:
>>>>>
>>>>> Aren't you just reinventing the inverse of a function?
>>>>>
>>>>
>>>> ---------------------------------------------------------------------------
>>>>>
>>>>> Jeff Newmiller The ..... ..... Go
>>>>
>>>> Live...
>>>>>
>>>>> DCN:<jdnewmil at dcn.davis.ca.us> Basics: ##.#. ##.#. Live
>>>>
>>>> Go...
>>>>>
>>>>> Live: OO#.. Dead: OO#..
>>>>
>>>> Playing
>>>>>
>>>>> Research Engineer (Solar/Batteries O.O#. #.O#. with
>>>>> /Software/Embedded Controllers) .OO#. .OO#.
>>>>
>>>> rocks...1k
>>>>>
>>>>>
>>>>
>>>> ---------------------------------------------------------------------------
>>>>>
>>>>> Sent from my phone. Please excuse my brevity.
>>>>>
>>>>> Emmanuel Levy <emmanuel.levy at gmail.com> wrote:
>>>>>
>>>>>> Hi,
>>>>>>
>>>>>> I am trying to normalize some data. First I fitted a principal curve
>>>>>> (using the LCPM package), but now I would like to apply a
>>>>>> transformation so that the curve becomes a "straight diagonal line"
>>>>
>>>> on
>>>>>>
>>>>>> the plot. The data used to fit the curve would then be normalized by
>>>>>> applying the same transformation to it.
>>>>>>
>>>>>> A simple solution could be to apply translations only (e.g., as done
>>>>>> after a fit using loess), but here rotations would have to be applied
>>>>>> as well. One could visualize this as the "stretching of a curve",
>>>>>> i.e., during such an "unfolding" process, both translations and
>>>>>> rotations would need to be applied.
>>>>>>
>>>>>> Before I embark on this (which would require me remembering long
>>>>>> forgotten geometry principles) I was wondering if you can think of
>>>>>> packages or tricks that could make my life simpler?
>>>>>>
>>>>>> Thanks for any input,
>>>>>>
>>>>>> Emmanuel
>>>>>>
>>>>>>
>>>>>> Below I provide an example - the black curve is to be "brought" along
>>>>>> the diagonal, and all data points normal to a "small segment" (of the
>>>>>> black curve) would undergo the same transformation as it - what
>>>>>> "small" is remains to be defined though.
>>>>>>
>>>>>> tmp=rnorm(2000)
>>>>>> X.1 = 5+tmp
>>>>>> Y.1 = 5+ (5*tmp+rnorm(2000))
>>>>>> tmp=rnorm(1000)
>>>>>> X.2 = 9+tmp
>>>>>> Y.2 = 40+ (1.5*tmp+rnorm(1000))
>>>>>> X.3 = 7+ 0.5*runif(500)
>>>>>> Y.3 = 15+20*runif(500)
>>>>>> X = c(X.1,X.2,X.3)
>>>>>> Y = c(Y.1,Y.2,Y.3)
>>>>>>
>>>>>> lpc1 = lpc(cbind(X,Y), scaled=FALSE, h=c(1,1) )
>>>>>> plot(lpc1)
>>>>>>
>>>>>> ______________________________________________
>>>>>> R-help at r-project.org mailing list
>>>>>> https://stat.ethz.ch/mailman/listinfo/r-help
>>>>>> PLEASE do read the posting guide
>>>>>> http://www.R-project.org/posting-guide.html
>>>>>> and provide commented, minimal, self-contained, reproducible code.
>>>>>
>>>>>
>>>
>>
>> ______________________________________________
>> R-help at r-project.org mailing list
>> https://stat.ethz.ch/mailman/listinfo/r-help
>> PLEASE do read the posting guide
>> http://www.R-project.org/posting-guide.html
>> and provide commented, minimal, self-contained, reproducible code.
>
>
> David Winsemius, MD
> West Hartford, CT
>
More information about the R-help
mailing list