aboutsummaryrefslogtreecommitdiff
path: root/pkg/images/immatch/src/geometry/geofunc.gx
diff options
context:
space:
mode:
Diffstat (limited to 'pkg/images/immatch/src/geometry/geofunc.gx')
-rw-r--r--pkg/images/immatch/src/geometry/geofunc.gx250
1 files changed, 250 insertions, 0 deletions
diff --git a/pkg/images/immatch/src/geometry/geofunc.gx b/pkg/images/immatch/src/geometry/geofunc.gx
new file mode 100644
index 00000000..3b34a207
--- /dev/null
+++ b/pkg/images/immatch/src/geometry/geofunc.gx
@@ -0,0 +1,250 @@
+include <math.h>
+include <math/gsurfit.h>
+
+$for (rd)
+
+# GEO_DROTMAG -- Adjust the coefficients of the fit using the database file.
+
+procedure geo_drotmag$t (dt, rec, sx1, sy1, xmag, ymag, xrot, yrot)
+
+pointer dt #I pointer to the text database file
+int rec #I record number
+pointer sx1, sy1 #I/O pointers to the x and y linear surfaces
+PIXEL xmag, ymag #I/O the x and y magnification
+PIXEL xrot, yrot #I/O the x and y axis rotation
+
+real dtgetr()
+
+begin
+ if (IS_$INDEF$T(xmag))
+ xmag = PIXEL (dtgetr (dt, rec, "xmag"))
+ if (IS_$INDEF$T(ymag))
+ ymag = PIXEL (dtgetr (dt, rec, "ymag"))
+ if (IS_$INDEF$T(xrot))
+ xrot = DEGTORAD (PIXEL(dtgetr (dt, rec, "xrotation")))
+ else
+ xrot = DEGTORAD(xrot)
+ if (IS_$INDEF$T(yrot))
+ yrot = DEGTORAD (PIXEL (dtgetr (dt, rec, "yrotation")))
+ else
+ yrot = DEGTORAD(yrot)
+ call geo_rotmag$t (sx1, sy1, xmag, ymag, xrot, yrot)
+end
+
+
+# GEO_DXYSHIFT -- Adjust the fitted xy shift using the database file.
+
+procedure geo_dxyshift$t (dt, rec, sx1, sy1, xout, yout, xref, yref,
+ xshift, yshift)
+
+pointer dt #I pointer to the text file database
+int rec #I the database record
+pointer sx1, sy1 #I/O pointers to the x and y linear surfaces
+PIXEL xout, yout #I the input coordinate system origin
+PIXEL xref, yref #I the reference coordinate system origin
+PIXEL xshift, yshift #I the origin shift in input coordinates
+
+$if (datatype == r)
+PIXEL gsgetr(), gseval()
+$else
+PIXEL dgsgetd(), dgseval()
+$endif
+
+begin
+$if (datatype == r)
+ if (IS_$INDEF$T(xref))
+ xref = (gsgetr (sx1, GSXMIN) + gsgetr (sx1, GSXMAX)) / 2.0
+ if (IS_$INDEF$T(yref))
+ yref = (gsgetr (sy1, GSYMIN) + gsgetr (sy1, GSYMAX)) / 2.0
+
+ if (IS_$INDEF$T(xout))
+ xout = gseval (sx1, xref, yref)
+ if (IS_$INDEF$T(yout))
+ yout = gseval (sy1, xref, yref)
+
+ if (IS_$INDEF$T(xshift))
+ xshift = xout - gseval (sx1, xref, yref)
+ if (IS_$INDEF$T(yshift))
+ yshift = yout - gseval (sy1, xref, yref)
+$else
+ if (IS_$INDEF$T(xref))
+ xref = (dgsgetd (sx1, GSXMIN) + dgsgetd (sx1, GSXMAX)) / 2.0d0
+ if (IS_$INDEF$T(yref))
+ yref = (dgsgetd (sy1, GSYMIN) + dgsgetd (sy1, GSYMAX)) / 2.0d0
+
+ if (IS_$INDEF$T(xout))
+ xout = dgseval (sx1, xref, yref)
+ if (IS_$INDEF$T(yout))
+ yout = dgseval (sy1, xref, yref)
+
+ if (IS_$INDEF$T(xshift))
+ xshift = xout - dgseval (sx1, xref, yref)
+ if (IS_$INDEF$T(yshift))
+ yshift = yout - dgseval (sy1, xref, yref)
+$endif
+
+ call geo_xyshift$t (sx1, sy1, xshift, yshift)
+end
+
+
+# GEO_ROTMAG -- Edit the coefficients of the linear surface which determine
+# magnification and rotation.
+
+procedure geo_rotmag$t (sx1, sy1, xscale, yscale, xrotation, yrotation)
+
+pointer sx1, sy1 #I/O pointers to the linear x and y surfaces
+PIXEL xscale, yscale #I the x and y scales
+PIXEL xrotation,yrotation #I the x and y axis rotation angles in radians
+
+PIXEL cosx, sinx, cosy, siny, xrange, yrange
+int ncoeff
+pointer sp, xcoeff, ycoeff
+$if (datatype == r)
+real gsgetr()
+int gsgeti()
+$else
+double dgsgetd()
+int dgsgeti()
+$endif
+
+begin
+ # Get the current solution.
+ call smark (sp)
+$if (datatype == r)
+ ncoeff = max (gsgeti (sx1, GSNSAVE), gsgeti (sy1, GSNSAVE))
+$else
+ ncoeff = max (dgsgeti (sx1, GSNSAVE), dgsgeti (sy1, GSNSAVE))
+$endif
+ call salloc (xcoeff, ncoeff, TY_PIXEL)
+ call salloc (ycoeff, ncoeff, TY_PIXEL)
+$if (datatype == r)
+ call gssave (sx1, Mem$t[xcoeff])
+ call gssave (sy1, Mem$t[ycoeff])
+$else
+ call dgssave (sx1, Mem$t[xcoeff])
+ call dgssave (sy1, Mem$t[ycoeff])
+$endif
+
+ # Define the scaling parameters.
+ cosx = cos (xrotation)
+ sinx = sin (xrotation)
+ cosy = cos (yrotation)
+ siny = sin (yrotation)
+
+ # Calculate coefficients.
+ Mem$t[xcoeff+GS_SAVECOEFF+1] = xscale * cosx
+ Mem$t[xcoeff+GS_SAVECOEFF+2] = yscale * siny
+ Mem$t[ycoeff+GS_SAVECOEFF+1] = -xscale * sinx
+ Mem$t[ycoeff+GS_SAVECOEFF+2] = yscale * cosy
+
+ # Normalize coefficients for-non polynomial functions.
+$if (datatype == r)
+ if (gsgeti (sx1, GSTYPE) != GS_POLYNOMIAL) {
+ xrange = gsget$t (sx1, GSXMAX) - gsget$t (sx1, GSXMIN)
+$else
+ if (dgsgeti (sx1, GSTYPE) != GS_POLYNOMIAL) {
+ xrange = dgsget$t (sx1, GSXMAX) - dgsget$t (sx1, GSXMIN)
+$endif
+ Mem$t[xcoeff+GS_SAVECOEFF+1] = Mem$t[xcoeff+GS_SAVECOEFF+1] *
+ xrange / 2.d0
+ Mem$t[xcoeff+GS_SAVECOEFF+2] = Mem$t[xcoeff+GS_SAVECOEFF+2] *
+ yrange / 2.d0
+ }
+$if (datatype == r)
+ if (gsgeti (sy1, GSTYPE) != GS_POLYNOMIAL) {
+ yrange = gsget$t (sy1, GSYMAX) - gsget$t (sy1, GSYMIN)
+$else
+ if (dgsgeti (sy1, GSTYPE) != GS_POLYNOMIAL) {
+ yrange = dgsget$t (sy1, GSYMAX) - dgsget$t (sy1, GSYMIN)
+$endif
+ Mem$t[ycoeff+GS_SAVECOEFF+1] = Mem$t[ycoeff+GS_SAVECOEFF+1] *
+ xrange / 2.d0
+ Mem$t[ycoeff+GS_SAVECOEFF+2] = Mem$t[ycoeff+GS_SAVECOEFF+2] *
+ yrange / 2.d0
+ }
+
+$if (datatype == r)
+ # Free the original fit.
+ call gsfree (sx1)
+ call gsfree (sy1)
+
+ # Restore the edited fit.
+ call gsrestore (sx1, Mem$t[xcoeff])
+ call gsrestore (sy1, Mem$t[ycoeff])
+$else
+ # Free the original fit.
+ call dgsfree (sx1)
+ call dgsfree (sy1)
+
+ # Restore the edited fit.
+ call dgsrestore (sx1, Mem$t[xcoeff])
+ call dgsrestore (sy1, Mem$t[ycoeff])
+$endif
+
+ call sfree (sp)
+end
+
+
+# GEO_XYSHIFT -- Shift the linear part of the fit in x and y.
+
+procedure geo_xyshift$t (sx1, sy1, xshift, yshift)
+
+pointer sx1, sy1 #I pointers to linear x and y surfaces
+PIXEL xshift, yshift #I the input x and y shifts
+
+int ncoeff
+pointer sp, xcoeff, ycoeff
+$if (datatype == r)
+int gsgeti()
+$else
+int dgsgeti()
+$endif
+
+begin
+ call smark (sp)
+
+ # Allocate working space.
+$if (datatype == r)
+ ncoeff = max (gsgeti (sx1, GSNSAVE), gsgeti (sy1, GSNSAVE))
+$else
+ ncoeff = max (dgsgeti (sx1, GSNSAVE), dgsgeti (sy1, GSNSAVE))
+$endif
+ call salloc (xcoeff, ncoeff, TY_PIXEL)
+ call salloc (ycoeff, ncoeff, TY_PIXEL)
+
+ # Get coefficients.
+$if (datatype == r)
+ call gssave (sx1, Mem$t[xcoeff])
+ call gssave (sy1, Mem$t[ycoeff])
+$else
+ call dgssave (sx1, Mem$t[xcoeff])
+ call dgssave (sy1, Mem$t[ycoeff])
+$endif
+
+ # Shift the coefficients.
+ Mem$t[xcoeff+GS_SAVECOEFF] = Mem$t[xcoeff+GS_SAVECOEFF] + xshift
+ Mem$t[ycoeff+GS_SAVECOEFF] = Mem$t[ycoeff+GS_SAVECOEFF] + yshift
+
+$if (datatype == r)
+ # Free original fit.
+ call gsfree (sx1)
+ call gsfree (sy1)
+
+ # Restore fit.
+ call gsrestore (sx1, Mem$t[xcoeff])
+ call gsrestore (sy1, Mem$t[ycoeff])
+$else
+ # Free original fit.
+ call dgsfree (sx1)
+ call dgsfree (sy1)
+
+ # Restore fit.
+ call dgsrestore (sx1, Mem$t[xcoeff])
+ call dgsrestore (sy1, Mem$t[ycoeff])
+$endif
+
+ call sfree (sp)
+end
+
+
+$endfor