Hmm, possibly the attached routine is of some help. It is extracted from a
thesis http://www.isical.ac.in/~amitava/thesis.ps
Regards
Karsten
Am Sat, 04 Feb 2006 21:28:07 +0100 schrieb Jeff N. <jnettle1@utk.edu>:
> Hi folks,
>
> I have a set of roughly circular ROI's for which I want to compute the
> maximum diameter - the longest straight line from one point on the
> boundary of the ROI, through the center, to the opposing boundary. Any
> ideas about the best way to do this? I do a loop I guess....compute
> the diameter for all 360 degrees around the center. Principal
> components analysis also springs to mind, but some of these ROI's have
> different scales in the x and y dimensions. Can anyone help?
>
> Jeff
>
; $Id: ad_msc.pro,v 1.2 2006/02/07 08:04:03 karo Exp $
;+
; NAME:
; AD_MSC
;
; PURPOSE:
; Calculation of the minimumspanning circle
;
; CATEGORY:
; Geometry
;
; CALLING SEQUENCE:
; Rad=ad_msc(X,Alpha)
;
; INPUTS:
; X: 2/3 dim Vector with coordinates [3,nx]
; Alpha: Weight
;
; KEYWORD PARAMETERS:
; WEIGHT: Switch set decreased weight, def. constant weight
;
; MAX_ITERATIONS: Number of iterations
; CENTER: variable, will contain the coordinates of msc
; EPS: Epsilon for exact match
; NP: Number of nonprinted lines if PRINT is on
; PRINT: Switch for printing
;
; OUTPUTS:
; Result: Radius of MSC
;
; PROCEDURE:
;
; EXAMPLE:
; See in program file
;
;
; MODIFICATION HISTORY:
; Written by: Karsten Rodenacker, after A. Datta, 1998.
; Nov. 99 Revised for 3d
; $Log: ad_msc.pro,v $
; Revision 1.2 2006/02/07 08:04:03 karo
; Conversion Win to Mac
;
;
FUNCTION ad_msc, x, a, $
MAX_ITERATIONS=nit, WEIGHT=ia, $
CENTER=cent,EPS=eps, NP=np, $
PRINT=print,HIT_POINT=i1
;-
IF NOT keyword_set(eps) THEN eps=0.1E-7
eps=double(eps)
; if not keyword_set(nit) then nit=2000l
IF NOT keyword_set(nit) THEN nit=(size(x,/dim))[1]*10.
IF NOT keyword_set(np) THEN np=100
sx=size(x)
nx=n_elements(x)/sx[1]
w=double(total(x[0,*])/nx)
FOR i=1,sx[1]-1 DO w=[temporary(w),double(total(x[i,*])/nx)]
it=0l
CASE sx[1] OF
2: REPEAT BEGIN
IF NOT keyword_set(ia) THEN alpha=double(a) $
ELSE alpha=double(a/(it+1))
mxd=max(sqrt((x[0,*]-w[0])^2+(x[1,*]-w[1])^2), i1)
pw=w
w=w+alpha*(x[*,i1]-w)
it=it+1
IF keyword_set(print) AND (it MOD np) THEN $
print,eps, ABS(pw-w)
ENDREP UNTIL (max(abs(pw-w)) LT eps) OR (it GE nit)
3: REPEAT BEGIN
IF NOT keyword_set(ia) THEN alpha=double(a) $
ELSE alpha=double(a/(it+1))
mxd=max(sqrt((x[0,*]-w[0])^2+(x[1,*]-w[1])^2+(x[2,*]-w[2])^2 ),
i1)
pw=w
w=w+alpha*(x[*,i1]-w)
it=it+1
IF keyword_set(print) AND (it MOD np) THEN $
print,eps, ABS(pw-w)
ENDREP UNTIL (max(abs(pw-w)) LT eps) OR (it GE nit)
ELSE: message,'dimension not implemented'
ENDCASE
cent=w
; return, sqrt((x[0,i1]-w[0])^2+(x[1,i1]-w[1])^2)
return, sqrt(total((x[*,i1]-w)^2))
END
; test wird nur bei .run ad_msc aufgerufen
x=randomn(seed,2,100)
rad=ad_msc(x,0.0001d,cent=c)
ci=ro_circle(c[0],c[1],rad)
window,0
plot,ci[0,*],ci[1,*]
plots,c[0],c[1],psym=2
oplot,x[0,*],x[1,*],/psym
x=randomn(seed,3,100)
rad=ad_msc(x,0.0001d,cent=c)
ci=ro_circle(c[0],c[1],rad)
window,1
plot_3dbox,x[0,*],x[1,*],x[2,*],psym=4,/solid
plots,ci[0,*],ci[1,*],replicate(c[2],32),col=0,/t3d
ci=ro_circle(c[1],c[2],rad)
plots,replicate(c[0],32),ci[0,*],ci[1,*],col=0,/t3d
ci=ro_circle(c[0],c[2],rad)
plots,ci[0,*],replicate(c[1],32),ci[1,*],col=0,/t3d
plots,[c[0]-rad,c[0]+rad],[c[1],c[1]],[c[2],c[2]],col=0,/t3d
plots,[c[0],c[0]],[c[1]-rad,c[1]+rad],[c[2],c[2]],col=0,/t3d
plots,[c[0],c[0]],[c[1],c[1]],[c[2]-rad,c[2]+rad],col=0,/t3d
wshow
!p.multi=0
END
--
Erstellt mit Operas revolutionï¿œrem E-Mail-Modul: http://www.opera.com/m2/
|
comp.lang.idl-pvwave archive
Members
Search
Help
Login
Home
