/*  このプログラムは，清水邦夫編著「地球環境データ」共立出版 (2002), pp.52-55
  に対応しています。このプログラムの利用条件は，GNU LGPL 

      http://www.gnu.org/licenses/licenses.html

  に準じます。

      This program is free software; you can redistribute it and/or modify
      it under the terms of the GNU General Public License as published by
      the Free Software Foundation; either version 2 of the License, or
      (at your option) any later version.

      This program is distributed in the hope that it will be useful,
      but WITHOUT ANY WARRANTY; without even the implied warranty of
      MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
      GNU General Public License for more details.

  Copyright 2002,2003 Shojiro Tanaka, Kenji Tahara                            */


#include <stdio.h>
#include <math.h>

#define   xyMAX   600           /* 数値地図の大きさ（近隣区画併合済データ，
                                   　ここでは正方形）  */
#define   uvMAX   1600          /* Landsat画像の大きさ */
#define   pixel   0.0285        /* ピクセルサイズ (km) */
#define   R       6377397.155   /* 地球の半径 (km)     */
#define   Z       705.0         /* 衛星高度 (km)       */
#define   pi      3.1415926535

int      x, y;   /* 数値地図座標 */
double   u, v;   /* 数値地図座標に対応する画像座標 */

int      m, p; 
double   n, q;

int      g1,g2,g3,g4; /* 周囲４隅Landsat画素値 */
double   gb;          /* 正射投影補正済画像値 */

double   mid,h,phi,SQsq,thetaS,gamma,phiS,deltaP;

/*精密補正用係数（GCPを用いて別途求める）*/
const double a1=1.999999,
             a2=2.999999,
             a3=3.999999,
             b1=4.999999,
             b2=5.999999,
             b3=6.999999;

FILE *fpR1, *fpR2, *fpW;

main()
{
   fpR1 = fopen("(Landsat ファイル名)", "rb"); /* Landsat原画像 */
   fpR2 = fopen("(数値地図ファイル名)", "r" ); /* 数値地図データ */
   fpW  = fopen("(補正出力ファイル名)", "w+"); /* 正射投影補正出力データ */
	
   u = 0.0;
   v = 0.0;

   for (y = 1; y <= xyMAX; y++)
   {
      for(x = 1; x <= xyMAX; x++)
      {
         /* 数値地図座標に対応する画像座標点の発生 */
         u = a1+a2*x+a3*y;
         v = b1+b2*x+b3*y;

         /* 衛星直下点のピクセル番号 */
         mid = 161.5*(v-77886.6)/2981.0;
			
         /* 標高を読取る */
         fscanf(fpR2,"%6d",&h);
         /* 海のデータ“-9999”を0にする */
         if (h < 0)
         {
            h = 0;
         }
         /* km単位系に変換する */
         h = h / 10000.0;

         /* 真の位置 */
         phi    = (mid-u)*pixel/R;
         SQsq   = (R+h)*(R+h)+(R+Z)*(R+Z)-2.0*(R+h)*(R+Z)*cos(phi);
         thetaS = asin( (R+h)*sin(phi)/sqrt(SQsq) );
         gamma  = pi - asin( (R+Z)*sin(thetaS)/R );
         phiS   = pi - thetaS - gamma;
 
        /* 起伏による位置ずれ（ピクセル）*/
         deltaP = R * (phiS - phi)/pixel;
         u      = u - deltaP;
         
         if (((u>=1) && (v>=1)) && ((uvMAX>=u) && (uvMAX>=v)))
         {
            uu = (int) u;
            vv = (int) v;
            fseek ( fpR1,uvMax*vv+uu, SEEK_SET );
            fscanf( fpR1,"%c",        &g1);
            fscanf( fpR1,"%c",        &g2);
            fseek ( fpR1,uvMAX-2,     SEEK_CUR );
            fscanf( fpR1,"%c",        &g3);
            fscanf( fpR1,"%c",        &g4);

            n = u - uu;
            q = v - vv;

            /* 共一次内挿，差が微少な場合，(255/2)*( (g_i-\mu_i)/\sigma_i + 1 ) 
               などどする規準化が必要な場合アリ */

            gb = ( g1*(1.0-r)*(1.0-s) + g2*r*(1.0-s) + g3*(1.0-r)*s + g4*r*s );
            if (gb > 255.0)
            {
               gb = 255.0;
            }
            else if (gb < 0.0)
            {
               gb =   0.0;
            }
            fprintf(fpW,"%4d",(int)gb);
            }

         /* 再配列できない */
         else
         {
            printf("ERROR");
         }
      }
      /* ラインの終わりに改行文字を入れる */
      fprintf(fpW,"%c",'\n');
   }
   fclose(fpR1);
   fclose(fpR2);
   fclose(fpW );

   return EXIT_SUCCESS;
} /* main */