Punti in quadrature su curve

Sergio_Alessi · 29 Gennaio 2022, 3:02pm

Mi dispiace Salvio ma non ho del codice condivisibile.

Sergio_Alessi · 29 Gennaio 2022, 3:36pm

Ciao Emilio!

certo che si riesce. Converti il polinomio in forma canonica,
calcoli derivata del polinomio e trovi gli zeri.

Sergio

0904 · 29 Gennaio 2022, 4:20pm

RiccardoMajewski:

Il secondo, esclude i casi dove manca la perpendicolarità tra tangente e vettore fornito:

  private void RunScript(Curve C, Vector3d V, ref object t)
  {
    double[] T = C.ExtremeParameters(V);
    List<double> T2 = new List<double>();
    foreach(double d in T){
      double a = Rhino.Geometry.Vector3d.VectorAngle(C.TangentAt(d), V);
      if(Math.Abs(a - Math.PI / 2) < this.RhinoDocument.ModelAbsoluteTolerance){
        T2.Add(d);
      }
    }
    t = T2;
  }

ma quindi, questo tipo di problema non rientra nella soluzione nell’uso del codice sopra citato?

RiccardoMajewski · 29 Gennaio 2022, 5:25pm

I “trucchi da strada” che faccio io sono accozzaglie di codice che imitano le cose imparate alle superiori… ma di fatto è tutto un po’ inventato alla fai-da-te …

Mi piacerebbe molto vedere un esempio minimale di quello che hai appena descritto…

Sergio_Alessi · 29 Gennaio 2022, 9:12pm

… però le applichi bene.

Curva bezier grado 3 con i seguenti punti di controllo

P0=(0.0,0.0)
P1=(1.0,1.0)
P2=(2.0,-1.0)
P3=(3.0;0.0)

'La curva ha equazione:

C(x,y)=P0*(1-t)^3+3*P1*t*(1-t)^2+3*P2*t^2*(1-t)+P3*t^3

dC/dt = 3*[(P1-P0)*t^2+2*(P2-P1)*t*(1-t)+(P3-P2)*t^2]

' si considera la componente Y (la X è monotona crescente ...)

dCy/dt = 3*[(1)*t^2+2*(-2)*t*(1-t)+(1)*t^2]

dCy/dt = 0

'tolgo il 3 perché è una costante ... risulta:

6*t^2-6*t+1=0

'si applica la nota formula ...

t=[6 +/- (36-24)^0.5]/(2*6)

luca.filippone · 30 Gennaio 2022, 3:28pm

Così troviamo gli zeri, ovvero i punti con derivata nulla rispetto al sistema di riferimento di Rhino, la funzione che ha mostrato Ricc li trova rispetto a qualsiasi sistema di riferimento (lo chiamava vettore perpendicolare alla tangente).

Sergio_Alessi · 30 Gennaio 2022, 4:03pm

Ciao Luca
questa è la soluzione per trovare i punti da snap quad.
Esiste la soluzione analitica anche per la condizione di perpendicolarità
punto - curva … ma anche curva - curva … anche tra superfici!
Le soluzioni sono analoghe.
Questa per me è materia studiata poco meno di trenta anni fa …
Lascio ai giovani come te di analisi I e II formulare la soluzione analitica.
Per concludere, non è detto che risolvere analiticamente il problema
sia la soluzione migliore. La soluzione brute di Emilio non è una cavolata, è di fatto
la più usata per molti problemi!

0904 · 30 Gennaio 2022, 11:46pm

prima di postare una domanda, faccio sempre qualche ricerca, ed in questo caso sono riuscito a trovare poco e nulla, credevo che sbagliassi le parole chiavi, ma a come vedo la cosa è un bel pò complessa.

RiccardoMajewski · 10 Febbraio 2022, 9:04pm

Sono riuscito a calcolare il punto di quadrante con il tuo metodo, matematicamente, con il metodo analitico come dicevi tu.

Ma non riesco comunque a capire questo passaggio:

Vedo che c’è la struttura del cubo di binomio (dove A=1-t e B=t):
A^3 + 3(A^2)B + 3A(B^2) + B^3
E la derivata:
3A^2 + 6AB + 3B^2

Ma non capisco il passaggio dai quattro punti P0,P1,P2,P3 ai 3 vettori(?) P1-P0,P2-P1,P3-P2 …
Però… funziona! … quindi ok.

Ma, potresti spiegarmi:

qual’è la funzione di una nurbs a grado 3 con 5+ punti di controllo?
pesi e knot? Immagino servirebbe integrarli per fare un lavoro completo…

Sergio_Alessi · 10 Febbraio 2022, 9:46pm

La derivata di una bezier si calcola cosi:
grado che moltiplica la bezier di grado n -1 con punti di controllo Pi - Pi-1.

Se hai una spline la decomponi in bezier con l’algoritmo knot insertion.
Se poi è razionale calcoli la derivata.

C = N / D
dC = (N’D -ND’)/D^2

emilio · 10 Febbraio 2022, 9:52pm

RiccardoMajewski:

Ma non riesco comunque a capire questo passaggio:
Sergio_Alessi:
> C(x,y)=P0*(1-t)^3+3*P1*t*(1-t)^2+3*P2*t^2*(1-t)+P3*t^3

dC/dt = 3*[(P1-P0)*t^2+2*(P2-P1)*t*(1-t)+(P3-P2)*t^2]
> ```

Grazie per l’info.

Ma (suppongo ovviamente) se derivi il polinomio e poi raccogli viene fuori la stessa cosa.

A me pero’ il primo termine viene in (1-t)^2 … giusto ?

Comunque derivando il polinomio il passaggio intermedio (se servisse) risulta …salvo errori :

dC/dt = 
-3*P0*(1-t)^2+
+3*P1*(1-t)^2+
-3*2*P1*t*(1-t)+
+3*2*P2*t*(1-t)+
-3*P2*t^2+
+3*P3*t^2

0904 · 10 Febbraio 2022, 10:02pm

un 3d aperto tempo fa proprio sulle Curve Nurbs,
bella tutta la discussione con tanti post interessanti.

Sergio_Alessi · 10 Febbraio 2022, 10:15pm

È giusto Emilio! Anche il passaggio intermedio.

Sergio_Alessi · 10 Febbraio 2022, 10:52pm

Condivido una soluzione in rhinoscript. Non l’ho testato molto. La precisione sembra buona
ma potrei sbagliarmi … Genera un errore con curve chiuse: vedrò con calma. Qui i limiti del
linguaggio iniziano a farsi sentire ciò nonostante non è lentissimo …
MinMaxCrv.rvb (5,4 KB)

Option Explicit
'Script written by Ing. Alessi Sergio
'Script copyrighted by Ing. Alessi Sergio
'Script version Alfa 0 date 02.10.2022

Call Main()
Sub Main()

	Dim crv:crv = rhino.getobject("Select a planar curve ...", 4)
	Dim g:g = rhino.CurveDegree(crv)
	Dim nK:nK = rhino.CurveKnotCount(crv)
	Dim K:K = rhino.CurveKnots(crv)
	Dim CPts():ReDim CPts((nk-g)*g*2)'Points coord.
	Dim PPts():ReDim PPts((nk-g)*g*2)'Points param.
	Dim Npts
	Dim d0, d1, i, j, h, par
	
	d0 = K(0)
	Cpts(0) = rhino.EvaluateCurve(crv, d0)
	Ppts(0) = d0
	h = 1
	i = g - 1
	While i < nk - 1
		d1 = K(i + 1)
		If d1 <> d0 Then
			For j=1 To 2 * g - 1
				par = d0 + (d1 - d0) * j / 2 / g
				PPts(h) = par
				CPts(h) = rhino.EvaluateCurve(crv, par)
				h = h + 1
			Next
			PPts(h) = d1
			CPts(h) = rhino.EvaluateCurve(crv, d1)
			d0 = d1
			h = h + 1
		End If
		i = i + 1
	Wend
	Npts = h
	
	Dim dd,ud, Iter,z
	Dim mMx():ReDim mMx(Npts)
	Dim mMy():ReDim mMy(Npts)
	Dim mME():ReDim mME(Npts)
	Dim ix
	Dim iy
	Dim aTol:aTol = 0.0001

	ix = 0
	iy = 0
	par = k(0)
	dd = rhino.CurveEvaluate(crv, par, 1)
	ud = rhino.VectorUnitize(dd(1))
	
	If abs(ud(0)) < aTol Then
		mMx(ix) = par
		ix = ix + 1
	End If
	
	If abs(ud(1)) < aTol Then
		mMy(iy) = par
		iy = iy + 1
	End If

	For i=1 To Npts - 2
		'********************************************************************************
		'MINMAX X
		'********************************************************************************
		If (CPts(i)(0) - CPts(i - 1)(0)) >= 0 And (CPts(i + 1)(0) - CPts(i)(0)) =< 0 Then
			par = PPts(i)
			dd = rhino.CurveEvaluate(crv, par, 1)
			ud = rhino.VectorUnitize(dd(1))
			If abs(ud(0)) > aTol Then
				If ud(0) > 0 Then
					d0 = PPts(i)
					d1 = PPts(i + 1)
				Else
					d0 = PPts(i - 1)
					d1 = PPts(i)
				End If
				Iter = 0
				While abs(ud(0)) > aTol And Iter < 50
					par = (d0 + d1) / 2
					dd = rhino.CurveEvaluate(crv, par, 1)
					ud = rhino.VectorUnitize(dd(1))
					If ud(0) > 0 Then
						d0 = Par
					Else
						d1 = Par
					End If
					iter = iter + 1
				Wend
				If iter >= 50 Then
					Call rhino.Print("Errore imprevisto param = " & par)
				Else
					mMx(ix) = par
					ix = ix + 1
				End If
			Else
				mMx(ix) = par
				ix = ix + 1
			End If
		End If
		
		If ((CPts(i)(0) - CPts(i - 1)(0)) < 0 And (CPts(i + 1)(0) - CPts(i)(0)) > 0) Then
			par = PPts(i)
			dd = rhino.CurveEvaluate(crv, par, 1)
			ud = rhino.VectorUnitize(dd(1))
			If abs(ud(0)) > aTol Then
				If ud(0) < 0 Then
					d0 = PPts(i)
					d1 = PPts(i + 1)
				Else
					d0 = PPts(i - 1)
					d1 = PPts(i)
				End If
				Iter = 0
				While abs(ud(0)) > aTol And Iter < 50
					par = (d0 + d1) / 2
					dd = rhino.CurveEvaluate(crv, par, 1)
					ud = rhino.VectorUnitize(dd(1))
					If ud(0) < 0 Then
						d0 = Par
					Else
						d1 = Par
					End If
					iter = iter + 1
				Wend
				If iter >= 50 Then
					Call rhino.Print("Errore imprevisto param = " & par)
				Else
					mMx(ix) = par
					ix = ix + 1
				End If
			Else
				mMx(ix) = par
				ix = ix + 1
			End If
		End If
		'********************************************************************************
		'MINMAX Y
		'********************************************************************************
		If (CPts(i)(1) - CPts(i - 1)(1)) >= 0 And (CPts(i + 1)(1) - CPts(i)(1)) =< 0 Then
			par = PPts(i)
			dd = rhino.CurveEvaluate(crv, par, 1)
			ud = rhino.VectorUnitize(dd(1))
			If abs(ud(1)) > aTol Then
				If ud(1) > 0 Then
					d0 = PPts(i)
					d1 = PPts(i + 1)
				Else
					d0 = PPts(i - 1)
					d1 = PPts(i)
				End If
				Iter = 0
				While abs(ud(1)) > aTol And Iter < 50
					par = (d0 + d1) / 2
					dd = rhino.CurveEvaluate(crv, par, 1)
					ud = rhino.VectorUnitize(dd(1))
					If ud(1) > 0 Then
						d0 = Par
					Else
						d1 = Par
					End If
					iter = iter + 1
				Wend
				If iter >= 50 Then
					Call rhino.Print("Errore imprevisto param = " & par)
				Else
					mMy(iy) = par
					iy = iy + 1
				End If
			Else
				mMy(iy) = par
				iy = iy + 1
			End If
		End If
		
		If ((CPts(i)(1) - CPts(i - 1)(1)) < 0 And (CPts(i + 1)(1) - CPts(i)(1)) > 0) Then
			par = PPts(i)
			dd = rhino.CurveEvaluate(crv, par, 1)
			ud = rhino.VectorUnitize(dd(1))
			If abs(ud(1)) > aTol Then
				If ud(1) < 0 Then
					d0 = PPts(i)
					d1 = PPts(i + 1)
				Else
					d0 = PPts(i - 1)
					d1 = PPts(i)
				End If
				Iter = 0
				While abs(ud(1)) > aTol And Iter < 50
					par = (d0 + d1) / 2
					dd = rhino.CurveEvaluate(crv, par, 1)
					ud = rhino.VectorUnitize(dd(1))
					If ud(1) < 0 Then
						d0 = Par
					Else
						d1 = Par
					End If
					iter = iter + 1
				Wend
				If iter >= 50 Then
					Call rhino.Print("Errore imprevisto param = " & par)
				Else
					mMy(iy) = par
					iy = iy + 1
				End If
			Else
				mMy(iy) = par
				iy = iy + 1
			End If
		End If

	Next

	par = k(nk - 1)
	dd = rhino.CurveEvaluate(crv, par, 1)
	ud = rhino.VectorUnitize(dd(1))

	If abs(ud(0)) < aTol Then
		mMx(ix) = par
		ix = ix + 1
	End If
	
	If abs(ud(1)) < aTol Then
		mMy(iy) = par
		iy = iy + 1
	End If
	
	rhino.EnableRedraw False
	For i=0 To  ix - 1
		par = rhino.AddPoint(rhino.EvaluateCurve(crv, mMx(i)))
	Next
	
	For i=0 To  iy - 1
		par = rhino.AddPoint(rhino.EvaluateCurve(crv, mMy(i)))
	Next
	rhino.EnableRedraw True
	
	
End Sub

Curva grado 7 con 275 punti di controllo.

emilio · 11 Febbraio 2022, 10:54am

Ciao Sergio

Sembra che il problema siano le curve periodiche.

Sergio_Alessi · 11 Febbraio 2022, 1:07pm

Vero Emilio

il vettore dei knot cambia … guarderò con molta calma.

Ciao!

emilio · 11 Febbraio 2022, 9:23pm

Scusa Sergio …

Se capisco come funziona lo script, all’inizio su ogni span della curva ricavi dei punti e ci costruisci una polilinea con 2*grado segmenti.
Poi dove la polilinea presenta un massimo o un minimo (rispetto alla coordinata scelta) cerchi il massimo o il minimo della curva.

Se e’ cosi’ … La mia curiosita’ e’ questa:
Come hai ricavato il numero 2*grado per i segmenti della polilinea ?

Sergio_Alessi · 12 Febbraio 2022, 12:20am

Ciao Emilio

che domande mi fai in pubblico!!! Il risultato è corretto?
Puoi tradurre il codice in py così mi fai sapere come gira? Con il poco tempo
a disposizione se lo scrivo in c++ ci metto un mese! … dimentico tutte le regole,
anche le basilari … . In VB scrivo di tutto …

emilio · 12 Febbraio 2022, 11:27am

Hahaha.

Ci provo .
( Se riesco a fare lo script lo posto cosi’ vedi tu come gira )
Come suggerisci di fare per le curve periodiche ?

Se e’ per GH puoi usare C#, per Rhino credo tocchi aspettare Rhino 8.

Sergio_Alessi · 12 Febbraio 2022, 2:04pm

Come faresti tu, vedo …
Il vettore dei knot cambia rispetto alle curve clamped.

Intendevo che il codice è palesemente orientato verso la programmazione con puntatori di memoria. C# come tutti i linguaggi net mmm …