1 #include <cmath>
2 #include <cstdlib>
3 #include <cstring>
4 #include <sstream>
5 using namespace std;
6 
7 #include "findFile.h"
8 #include "xpUtil.h"
9 
10 #include "libimage/Image.h"
11 
12 /*
13   Convert lon, lat coordinates to screen coordinates for the Mollweide
14   projection.
15 */
16 static void
sphericalToPixel(const int width,const int height,const double lon,const double lat,int & x,int & y)17 sphericalToPixel(const int width, const int height,
18                  const double lon, const double lat,
19                  int &x, int &y)
20 {
21     double theta = lat;
22     double del_theta = 1;
23     while (fabs(del_theta) > 1e-5)
24     {
25         del_theta = -((theta + sin(theta) - M_PI * sin(lat))
26                       / (1 + cos(theta)));
27         theta += del_theta;
28     }
29     theta /= 2;
30 
31     double X = lon / M_PI * cos(theta);
32     double Y = sin(theta);
33 
34     x = (int) ((X + 1) * width/2);
35     y = (int) (height/2 * (1 - Y));
36 }
37 
38 /*
39   Convert x, y coordinates on the screen to lon, lat for a rectangular
40   projection.
41 */
42 static void
pixelToSpherical(const int width,const int height,const int x,const int y,double & lon,double & lat)43 pixelToSpherical(const int width, const int height,
44                  const int x, const int y,
45                  double &lon, double &lat)
46 {
47     lon = (x + 0.5) * TWO_PI / width - M_PI;
48     lat = M_PI_2 - (y + 0.5) * M_PI / height;
49 }
50 
51 static void
equalizeHistogram(unsigned char * & rgb,const int width,const int height)52 equalizeHistogram(unsigned char *&rgb, const int width, const int height)
53 {
54     // create a histogram
55     int *hist = new int[256];
56     for (int i = 0; i < 256; i++) hist[i] = 0;
57     for (int i = 0; i < 3 * width * height; i += 3)
58         hist[(int) rgb[i]]++;
59 
60     // create an integrated histogram
61     int *ihist = new int[256];
62     ihist[0] = hist[0];
63     for (int i = 1; i < 256; i++) ihist[i] = ihist[i-1] + hist[i];
64 
65     // replace the histogram by an intensity map
66     double denom = static_cast<double> (ihist[255] - ihist[0]);
67     if (denom > 0)
68     {
69         for (int i = 0; i < 256; i++)
70             hist[i] = static_cast<int> (255 * (ihist[i] - ihist[0])/denom);
71 
72         for (int i = 0; i < 3 * width * height; i++)
73             rgb[i] = static_cast<unsigned char> (hist[static_cast<int>(rgb[i])]);
74     }
75 
76     delete [] hist;
77     delete [] ihist;
78 }
79 
80 /*
81   This routine reads in a global cloud image downloaded from
82   http://www.ssec.wisc.edu/data/comp/latest_moll.gif
83   and reprojects and resizes the image, gets rid of the ugly pink
84   coastlines, and stretches the contrast.
85 */
86 static bool
convertSsecImage(Image * image,unsigned char * & rgb)87 convertSsecImage(Image *image, unsigned char *&rgb)
88 {
89     // There's a 20 pixel border on the left & right and a 10 pixel border
90     // on the top and bottom
91     image->Crop(10, 20, image->Width() - 10, image->Height() - 20);
92     const int image_height = image->Height();
93     const int image_width = image->Width();
94 
95     // This array will hold the final cloud image
96     rgb = (unsigned char *) malloc(3 * image_width * image_height);
97     if (rgb == NULL) return(false);
98 
99     memset(rgb, 0, 3 * image_width * image_height);
100 
101     // This converts the mollweide projection to rectangular
102     double lon, lat;
103     int ipos = 0;
104     const unsigned char *moll = image->getRGBData();
105     for (int j = 0; j < image_height; j++)
106     {
107         for (int i = 0; i < image_width; i++)
108         {
109             int ii, jj;
110             pixelToSpherical(image_width, image_height, i, j, lon, lat);
111             sphericalToPixel(image_width, image_height, lon, lat, ii, jj);
112             memcpy(rgb + ipos, moll + 3 * (jj * image_width + ii), 3);
113             ipos += 3;
114         }
115     }
116 
117     int avg;
118     int npoints;
119     int avgwhole = 0;
120     int npointswhole = 0;
121 
122     // Replace pink outlines by the average value in a 10x10 pixel square.
123     for (int j = 0; j < 31; j++)
124     {
125         for (int i = 0; i < 62; i++)
126         {
127             avg = 0;
128             npoints = 0;
129             for (int jj = 0; jj < 10; jj++)
130             {
131                 for (int ii = 0; ii < 10; ii++)
132                 {
133                     ipos = 3*((10 * j + jj) * 620 + 10 * i + ii);
134                     if (!(rgb[ipos] == 0xff
135                           && rgb[ipos+1] == 0
136                           && rgb[ipos+2] == 0xff))
137                     {
138                         npoints++;
139                         avg += (int) rgb[ipos];
140                         npointswhole++;
141                         avgwhole += (int) rgb[ipos];
142                     }
143                 }
144             }
145             if (npoints != 0) avg = (int) (avg / (double) npoints);
146 
147             for (int jj = 0; jj < 10; jj++)
148             {
149                 for (int ii = 0; ii < 10; ii++)
150                 {
151                     ipos = 3*((10 * j + jj) * 620 + 10 * i + ii);
152                     if (rgb[ipos] == 0xff
153                         && rgb[ipos+1] == 0
154                         && rgb[ipos+2] == 0xff)
155                         memset(rgb + ipos, avg, 3);
156                 }
157             }
158         }
159     }
160     // Fill in the poles
161     if (npointswhole != 0)
162         avgwhole = (int) (avgwhole / (double) npointswhole);
163     for (int i = 0; i < image_width * image_height; i++)
164     {
165         ipos = 3 * i;
166         if (rgb[ipos] < 0x03) memset(rgb + ipos, avgwhole, 3);
167     }
168 
169     // Smooth out the seam at 180 degrees Longitude
170     double eastVal, westVal;
171     int eastIndex, westIndex;
172     for (int i = 0; i < image_height - 1; i++)
173     {
174         eastIndex = 3 * (i * image_width + 1);
175         westIndex = 3 * ((i + 1) * image_width - 2);
176         eastVal = (double) rgb[eastIndex];
177         westVal = (double) rgb[westIndex];
178         memset(rgb + eastIndex - 3,
179                (int) (eastVal + (westVal - eastVal)/3), 3);
180         memset(rgb + westIndex + 3,
181                (int) (westVal + (eastVal - westVal)/3), 3);
182     }
183 
184     equalizeHistogram(rgb, image_width, image_height);
185 
186     return(true);
187 }
188 
189 void
loadSSEC(Image * & image,const unsigned char * & rgb,string & imageFile,const int imageWidth,const int imageHeight)190 loadSSEC(Image *&image, const unsigned char *&rgb, string &imageFile,
191          const int imageWidth, const int imageHeight)
192 {
193     bool foundFile = findFile(imageFile, "images");
194     if (foundFile)
195     {
196         image = new Image;
197         foundFile = image->Read(imageFile.c_str());
198     }
199 
200     if (foundFile)
201     {
202         unsigned char *tmpRGB = NULL;
203         if (!convertSsecImage(image, tmpRGB))
204         {
205             ostringstream errStr;
206             errStr << "Can't read SSEC map file: " << imageFile << "\n";
207             xpWarn(errStr.str(), __FILE__, __LINE__);
208             return;
209         }
210 
211         Image *tmpImage = image;
212         image = new Image(image->Width(), image->Height(), tmpRGB, NULL);
213         delete tmpImage;
214         free(tmpRGB);
215         if ((image->Width() != imageWidth)
216             || (image->Height() != imageHeight))
217         {
218             ostringstream errStr;
219             errStr << "Resizing SSEC cloud map\n"
220                    << "For better performance, all image maps should "
221                    << "be the same size as the day map\n";
222             xpWarn(errStr.str(), __FILE__, __LINE__);
223             image->Resize(imageWidth, imageHeight);
224         }
225         rgb = image->getRGBData();
226     }
227     else
228     {
229         ostringstream errStr;
230         errStr << "Can't load map file " << imageFile << "\n";
231         xpWarn(errStr.str(), __FILE__, __LINE__);
232     }
233 }
234 
235