PROGRAM xdigit * Calculate x, y coordinates from digitizer readout (u, v) * ディジタイザ(Graphtec)のよみ (u, v) から (x, y)座標を計算する * 1995-10-03, 1996-12-16, 2000-06-29 by K. Masuda * アフィン変換 x = a u + b v + c * y = d u + e v + f * を、コントロールポイントの座標値を規格化したものに適用 * (規格化のため、c = 0, f = 0 となる)。 * IMPLICIT NONE INTEGER nctldm, nobsdm * 配列の寸法はデータが読み込めるようにとること PARAMETER (nctldm=50, nobsdm=1000) CHARACTER*40 hfctl, hfobs, hfout REAL xctl(nctldm),yctl(nctldm),uctl(nctldm),vctl(nctldm) REAL xobs(nobsdm),yobs(nobsdm),uobs(nobsdm),vobs(nobsdm) REAL xtmp,ytmp,utmp,vtmp,dummy,dobs(nobsdm) INTEGER nctl, kctl, nobs, kobs REAL*8 xmean, ymean, umean, vmean REAL*8 xsigma, ysigma, usigma, vsigma REAL*8 suv, sux, suy, svx, svy REAL*8 aa, bb, dd, ee * * BEGIN * ファイル名の決定 * CALL getarg(1, hfctl) CALL getarg(2, hfobs) CALL getarg(3, hfout) cc WRITE (*,*) 'file name for control points - ', hfctl cc WRITE (*,*) 'file name for observation points - ', hfobs cc WRITE (*,*) 'file name for output - ', hfout * * コントロールポイントの処理 * データを読む (x, y座標をprojで追加したもの) * OPEN(unit=11,file=hfctl,status='OLD', & access='SEQUENTIAL',form='FORMATTED') kctl = 0 9000 CONTINUE READ(11,*,END=9900) xtmp,ytmp,dummy,utmp,vtmp kctl = kctl+1 uctl(kctl) = utmp vctl(kctl) = vtmp xctl(kctl) = xtmp yctl(kctl) = ytmp GO TO 9000 9900 CONTINUE CLOSE(unit=11) nctl = kctl * * 平均 * xmean = 0.0d0 ymean = 0.0d0 umean = 0.0d0 vmean = 0.0d0 DO kctl = 1, nctl xmean = xmean + xctl(kctl) ymean = ymean + yctl(kctl) umean = umean + uctl(kctl) vmean = vmean + vctl(kctl) END DO xmean = xmean / nctl ymean = ymean / nctl umean = umean / nctl vmean = vmean / nctl * * 標準偏差 * xsigma = 0.0d0 ysigma = 0.0d0 usigma = 0.0d0 vsigma = 0.0d0 DO kctl = 1, nctl xsigma = xsigma + (xctl(kctl)-xmean)**2 ysigma = ysigma + (yctl(kctl)-ymean)**2 usigma = usigma + (uctl(kctl)-umean)**2 vsigma = vsigma + (vctl(kctl)-vmean)**2 END DO xsigma = sqrt( xsigma / nctl ) ysigma = sqrt( ysigma / nctl ) usigma = sqrt( usigma / nctl ) vsigma = sqrt( vsigma / nctl ) * ** WRITE (*,*) xmean, ymean, umean, vmean ** WRITE (*,*) xsigma, ysigma, usigma, vsigma * * データを規格化 * DO kctl = 1, nctl xctl(kctl) = (xctl(kctl) - xmean) / xsigma yctl(kctl) = (yctl(kctl) - ymean) / ysigma uctl(kctl) = (uctl(kctl) - umean) / usigma vctl(kctl) = (vctl(kctl) - vmean) / vsigma ** WRITE (*,*) xctl(kctl), yctl(kctl), ** & uctl(kctl), vctl(kctl) END DO * * 共分散 * suv = 0.0d0 sux = 0.0d0 suy = 0.0d0 svx = 0.0d0 svy = 0.0d0 DO kctl = 1, nctl suv = suv + uctl(kctl)*vctl(kctl) sux = sux + uctl(kctl)*xctl(kctl) suy = suy + uctl(kctl)*yctl(kctl) svx = svx + vctl(kctl)*xctl(kctl) svy = svy + vctl(kctl)*yctl(kctl) END DO suv = suv / nctl sux = sux / nctl suy = suy / nctl svx = svx / nctl svy = svy / nctl * * アフィン変換の係数を求める * aa = (sux - suv*svx) / (1.0d0 - suv**2) bb = (svx - suv*sux) / (1.0d0 - suv**2) dd = (suy - suv*svy) / (1.0d0 - suv**2) ee = (svy - suv*suy) / (1.0d0 - suv**2) ** WRITE (*,*) aa, bb ** WRITE (*,*) dd, ee * * 観測点の処理 * データを読む * OPEN(unit=11,file=hfobs,status='OLD', & access='SEQUENTIAL',form='FORMATTED') kobs = 0 19000 CONTINUE READ(11,*,END=19900) utmp,vtmp,dummy kobs = kobs+1 uobs(kobs) = utmp vobs(kobs) = vtmp dobs(kobs) = dummy GO TO 19000 19900 CONTINUE CLOSE(unit=11) nobs = kobs * * 換算と出力 * OPEN(unit=11,file=hfout,status='UNKNOWN', & access='SEQUENTIAL',form='FORMATTED') DO kobs = 1, nobs * * コントロールポイントの(u, v)の平均と標準偏差を使って規格化 * uobs(kobs) = (uobs(kobs) - umean) / usigma vobs(kobs) = (vobs(kobs) - vmean) / vsigma * * アフィン変換 * xobs(kobs) = aa * uobs(kobs) + bb * vobs(kobs) yobs(kobs) = dd * uobs(kobs) + ee * vobs(kobs) * * 規格化の逆算 * xobs(kobs) = xsigma * xobs(kobs) + xmean yobs(kobs) = ysigma * yobs(kobs) + ymean * * 結果を書き出す (フォーマットは内容によって変えること) * WRITE (11,'(3F12.1)') xobs(kobs), yobs(kobs), dobs(kobs) END DO CLOSE(unit=11) END