C PROGRAM BNDBE.FT C ---------------- C C SUBROUTINE BNDBE(BE,P,ICH1,ICH2,MDIFF,AVC,AVD,KS,IWIND) C C PETER LEMKIN C IMAGE PROCESSING UNIT, DCBD C NATIONAL CANCER INSTITUTE C NATIONAL INSTITUTES OF HEALTH C 9000 ROCKVILLE PIKE C BETHESDA, MD. 20014 C C JULY 14, 1977 C JULY 12, 1977 C JULY 11, 1977 C JULY 7, 1977 C JULY 5, 1977 C JULY 1, 1977 C JUNE 26, 1977 C JUNE 23, 1977 C C INTRODUCTION C ------------ C BNDBE COMPUTES THE BENDING ENERGY FROM CHAIN CODE C DIFFERENCES AS FOLLOWS: C C BE= SUM (1/PERIMETER)*(CURVATURE/UNITLENGTH)**2 C C ARGUMENTS C --------- C BE - COMPUTED BENDING ENERGY (ZERO EXTERNALLY) C P - TOTAL PERIMETER (SUPPLIED TO BNDBE) C ICH1 - CURRENT CHAIN CODE (SUPPLIED TO BNDBE) C ICH2 - PREVIOUE CHAIN CODE (SUPPLIED TO BNDBE) C MDIFF - COMPUTED CHAIN CODE DIFFERENCE (-4 TO +4) C AVC - COMPUTED AVERAGE CHAIN CODE FOR WIDTH IWIND C CENTERED AT POINT (KS-IWIND/2). C NOTE: AVC FOR POINT IWIND RETURNED! C AVD - COMPUTED CHAIN CODE DIFFERENCE CENTERED AT C POINT IWIND: C (A1-A2), C WHERE: A1=ATAN(IAY1/IAX1) C WHERE: IAY1=SUM IAYQ(i), IAX1=SUM IAXQ(i) C FOR i=KS+1+IWIND/2:KS+IWIND, C A2=ATAN(IAY2/IAX2) C WHERE: IAY2=SUM IAYQ(i), IAX2=SUM IAXQ(i) C FOR i=KS+IWIND/2+1:KS+IWIND(3/2). C NOTE: AVCi'<==(If (AVCi > 4.0) C Then AVCi C Else AVCi-8.0); C NOTE: AVD FOR POINT (KS+IWIND-1) RETURNED. C KS - CURRENT POINT NUMBER C IWIND - WINDOW SIZE C C RELATION OF CHAIN CODE TO AVC AND AVD C ------------------------------------- C C NOTE: FOR WINDOW [i:i+W] C The chain code queues ICHXQ, ICHYQ holds data from [i:i+W]. C C The ACQ (average chain code queue) holds data from C [i+W/2:i+(3/2)W]. C IAXQ, IAYQ hold the average x,y coordinates used C to form AVCi by ATAN function. C C The AVD (the 1/2 difference of the ACQ) holds at (i+(3/2)W). C OPDEFS C ------ S OPDEF TADI 1400 S OPDEF ISZI 2400 S OPDEF DCAI 3400 C C S OPDEF MQA 7501 S OPDEF MQL 7421 S OPDEF KRS 6034 S OPDEF BSW 7002 C C S OPDEF LDXP 6443 S OPDEF LDYP 6444 C S OPDEF DISP1 6435 S OPDEF DISP2 6436 C C DIMENSION ACQ(13) C DIMENSION IAXQ(13),IAYQ(13) C DIMENSION ICHXQ(13),ICHYQ(13) C [1] DEFINE CHAIN DIFFERENCE: ICDIFF KDIFF=ICH1-ICH2 C COMPUTE: ICDIFF=IABS(KDIFF) S TAD \KDIFF S SPA S CIA S DCA \ICDIFF C C IF (ICDIFF>4) C THEN ICDIFF<==8-ICDIFF, C IF KDIFF > 0 C THEN KDIFF<==8-KDIFF C ELSE KDIFF<==KDIFF+8; S TAD \ICDIFF S TAD (-5 S SPA CLA S JMP \102 /OK, CONTINUE C C THEN INVERT! ICDIFF=8-ICDIFF C S TAD \KDIFF S SPA CLA S JMP \103 /ELSE KDIFF=8-KDIFF GOTO 102 103 KDIFF=-(KDIFF+8) C PASS IT BACK BY REFERENCE 102 MDIFF=KDIFF C C C [2] DEFINE PATH LENGTH: FLENGTH C NOTE: PIXEL DISTANCES C FLENGTH <== (JC2_ICH2 LAND '1)*0.207+ C (JC1_ICH1 LAND '1)*0.207+1.0; S TAD I \ICH2 S AND (0001 S DCA \JC2 S TAD I \ICH1 S AND (0001 S DCA \JC1 FLENGTH=1.0+0.207*FLOAT(JC2+JC1) C C [3] DEFINE CURVATURE CURVATURE: CURVATURE=FLOAT(ICDIFF)/FLENGTH C C C [4] COMPUTE INCREMENTAL AND TOTAL BENDING ENERGY CURVATURE=(CURVATURE*CURVATURE)/P BE=BE+CURVATURE C C C [5] IF IWIND =1 C THEN RETURN S CLA CMA S TAD I \IWIND S SNA CLA S JMP \2047 C C C ELSE COMPUTE THE RUNNING AVERAGES. C IF KS=1 THEN FIRST POINT , INIT AVGS C C CVT CHAIN CODE (ICH1) TO (IX,IY) DIFFERENCE COORDINATES S JMS CVTCHXY C IF(KS-1)510,501,510 C C C [5.1] INITIALIZE FIRST PIXEL C ZERO QUEUES 501 CONTINUE C C COMPUTE THE SCALE FACTOR TO CVT + RADIANS TO CHAIN CODE C NOTE: PI RADIANS = 180 DEGREES C NOTE: CHAIN CODE = ANGLE/45 DEGREES. C NOTE: SC=180.0/(45.0*3.14159265) SC=1.273240 C C FLOAT THE SIZE KSIZE=MIN(IWIND,13) ISIZM1=KSIZE-1 C C COMPUTE: KSZHALF=(KSIZE-3)/2 S TAD (-3 S TAD \KSIZE S CLL RAR /DIVIDE BY 2 S DCA \KSZHALF C C C [5.2] ACCUMULATE DATA, REMEMBER N'TH DATA POINT C IN QUEUES C IF KSIZE > KS C THEN K<==KS C ELSE K<==KSIZE 510 K=KSIZE IF(KS-KSIZE)525,524,524 525 K=KS C C C C [5.2.1] ADVANCE ENTRY ONTO THE QUEUES C 524 IXOLD=0 IYOLD=0 IAX1=0 IAX2=0 IAY1=0 IAY2=0 C DO 511 I=1,ISIZM1 C J IS POINTER TO SINK , J-1 PTR TO SOURCE C THUS: J RANGES [KSIZE:2], J-1 TO [KSIZE-1:1]. J=KSIZE+1-I C C MOVE THE AVG CHAIN VALUE DOWN THE QUEUE C COMPUTE: AA=ACQ(J-1) JPTR=(J-1)-1 JPTR=JPTR+JPTR+JPTR S TAD \JPTR S TAD PACQ S DCA GA1# S TAD GA1# S TAD (3 S DCA GA2# S CALL 1,FAD S GA1, ARG \ACQ /WILL CONTAIN ACQ(J-1) S CALL 1,STO S ARG \AA C C COMPUTE: ACQ(J)=AA S CALL 1,FAD S ARG \AA S CALL 1,STO S GA2, ARG \ACQ /WILL CONTAIN ACQ(J); C C C C MOVE IAXQ(J)<==IAX<==IAXQ(J-1) C COMPUTE: IAX=IAXQ(J-1) S TAD (-3 S TAD \J S TAD PIAXQ S DCA 11 S TADI 11 S DCA \IAX C C COMPUTE: IAXQ(J)=IAX S TAD \IAX S DCAI 11 C C C MOVE IAYQ(J)<==IAY<==IAYQ(J-1) C COMPUTE: IAY=IAYQ(J-1) S TAD (-3 S TAD \J S TAD PIAYQ S DCA 11 S TADI 11 S DCA \IAY C C COMPUTE: IAYQ(J)=IAY S TAD \IAY S DCAI 11 C C C IF 1 LEQ J LEQ ((KSIZE-3)/2) C THEN IAX1_IAX1+IAX, IAY1_IAY1+IAY IF(J-KSZHALF)527,527,528 527 IAX1=IAX1+IAX IAY1=IAY1+IAY GOTO 529 C C IF ((KSIZE(KSIZE-3)/2)+4) LEQ J LEQ KSIZE C THEN IAX2_IAX2+IAX, IAY2_IAY2+IAY 528 IF(J-KSZHALF+4)529,526,526 526 IAX2=IAX2+IAX IAY2=IAY2+IAY 529 CONTINUE C C C C MOVE THE X STACK AND GET X VALUE C COMPUTE: IXX=ICHXQ(J-1) S TAD (-3 S TAD \J S TAD PICHXQ S DCA 11 S TADI 11 S DCA \IXX C C COMPUTE: ICHXQ(J)=IXX S TAD \IXX S DCAI 11 C C MOVE THE Y STACK AND GET Y VALUE C COMPUTE: IYY=ICHYQ(J-1) S TAD (-3 S TAD \J S TAD PICHYQ S DCA 11 S TADI 11 S DCA \IYY C C COMPUTE: ICHYQ(J)=IYY S TAD \IYY S DCAI 11 C C C IF J LEQ K S TAD \J S CIA S TAD \K S SPA CLA S JMP \511 /DON'T ADD IT C THEN ADD IXX, IYY IXOLD=IXOLD+IXX IYOLD=IYOLD+IYY 511 CONTINUE C C C [5.2.2] COMPUTE AVG (X,Y) VALUES AND QUEUE THEM IAXQ=(IXOLD+IX) IAYQ=(IYOLD+IY) C C C [5.2.3] PUSH NEW DATA ONTO FRONT OF QUEUE ICHXQ=IX ICHYQ=IY C C C [5.3] RETURN AVG CHAIN CODE AT ((KSIZE-3)/2)+2) AVC=ACQ(KSZHALF+2) C IF AVC < 0.0 C THEN AVC<==AVC+8.0; IF(AVC)531,530,530 531 AVC=AVC+8.0 C C C [5.3.1] NOTE: RETURN A DECIMAL CHAIN CODE DIRECTION C NOTE: ARCTAN RETURNS + RADIANS 530 AA=SC*ARCTAN(FLOAT(IAYQ),FLOAT(IAXQ)) C C PUSH INTO FRONT OF QUEUE C C COMPUTE: ACQ(1)=AA S CALL 1,FAD S ARG \AA S CALL 1,STO S ARG \ACQ C C UPDATE A1 IAX1=IAX1+IAXQ IAY1=IAY1+IAYQ C C C [5.4] COMPUTE AVERAGE DIFFERENCE: C COMPUTE THE TWO ANGLES A1=SC*ARCTAN(FLOAT(IAY1),FLOAT(IAX1)) A2=SC*ARCTAN(FLOAT(IAY2),FLOAT(IAX2)) C CALL BNDHEMISPHERE(A1,A2,ADIFF) C AVD=ADIFF C C C@C***DEBUG*** C@ WRITE(3,1234)KS,IX,IY,A1,A2,AVC,AVD C@1234 FORMAT(/,10X,'KS=',I3,'[',I3,':',I3,']',', A1=',F7.2,', A2=' C@ 1,F7.2,', AVC=',F7.2,', AVD=',F7.2) C@ DO 1235 I=1,K C@ IXX=ICHXQ(I) C@ IYY=ICHYQ(I) C@ AAA=ACQ(I) C@ IAX=IAXQ(I) C@ IAY=IAYQ(I) C@1235 WRITE(3,1236)I,IXX,IYY,AAA,IAX,IAY C@1236 FORMAT(10X,'#',I3,', XQ=',I5,', YQ=',I5,', ACQ=',F7.3 C@ 1,', IAX=',I5,', IAY=',I5) C@C************ 2047 RETURN C ******************************************************** C SUBROUTINE C V T C H X Y C ********************************************************* C CONVERT ICH1 TO IX,IY DIFFERENCE COORDINATES S CPAGE 3 S RCVTCHXY, JMP I CVTCHXY S CVTCHXY, 0 M=ICH1 S TAD \M S AND (0007 S IAC S DCA \M GOTO(1000,1001,1002,1003,1004,1005,1006,1007),M C 1000 IX=+1 IY=0 S JMP RCVTCHXY C 1001 IX=+1 IY=+1 S JMP RCVTCHXY C 1002 IX=0 IY=+1 S JMP RCVTCHXY C 1003 IX=-1 IY=+1 S JMP RCVTCHXY C 1004 IX=-1 IY=0 S JMP RCVTCHXY C 1005 IX=-1 IY=-1 S JMP RCVTCHXY C 1006 IX=0 IY=-1 S JMP RCVTCHXY C 1007 IX=+1 IY=-1 S JMP RCVTCHXY C C C ****POINTERS*** S PICHYQ, \ICHYQ S PICHXQ, \ICHXQ S PACQ, \ACQ S PIAXQ, \IAXQ S PIAYQ, \IAYQ S PAA, \AA END