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 <dwinsem...@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 <jdnew...@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:<jdnew...@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.l...@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 <jdnew...@dcn.davis.ca.us>
>>>> wrote:
>>>>>
>>>>> Aren't you just reinventing the inverse of a function?
>>>>>
>>>>
>>>> ---------------------------------------------------------------------------
>>>>>
>>>>> Jeff Newmiller The ..... ..... Go
>>>>
>>>> Live...
>>>>>
>>>>> DCN:<jdnew...@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.l...@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@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@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
>
______________________________________________
R-help@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.
