1 //*******************************************************************
2 // Copyright (C) 2000 ImageLinks Inc.
3 //
4 // License: LGPL
5 //
6 // See LICENSE.txt file in the top level directory for more details.
7 //
8 // Author: Garrett Potts
9 //
10 //*************************************************************************
11 // $Id: ossimFeatherMosaic.cpp 15766 2009-10-20 12:37:09Z gpotts $
12
13 #include <ossim/imaging/ossimFeatherMosaic.h>
14 #include <ossim/base/ossimDpt.h>
15 #include <ossim/imaging/ossimImageData.h>
16 #include <ossim/base/ossimConstants.h>
17 #include <ossim/base/ossimLine.h>
18 #include <ossim/base/ossimTrace.h>
19
20 static ossimTrace traceDebug("ossimFeatherMosaic:debug");
21
22 RTTI_DEF1(ossimFeatherMosaic, "ossimFeatherMosaic", ossimImageMosaic);
23
24 using namespace std;
25
ossimFeatherMosaic()26 ossimFeatherMosaic::ossimFeatherMosaic()
27 :ossimImageMosaic(),
28 theInputFeatherInformation(NULL),
29 theAlphaSum(NULL),
30 theResult(NULL),
31 theFeatherInfoSize(0)
32 {
33 }
34
ossimFeatherMosaic(ossimConnectableObject::ConnectableObjectList & inputSources)35 ossimFeatherMosaic::ossimFeatherMosaic(ossimConnectableObject::ConnectableObjectList& inputSources)
36 :ossimImageMosaic(inputSources),
37 theInputFeatherInformation(NULL),
38 theAlphaSum(NULL),
39 theResult(NULL),
40 theFeatherInfoSize(0)
41 {
42 initialize();
43 }
44
~ossimFeatherMosaic()45 ossimFeatherMosaic::~ossimFeatherMosaic()
46 {
47 if(theInputFeatherInformation)
48 {
49 delete [] theInputFeatherInformation;
50 theInputFeatherInformation = NULL;
51 }
52 theFeatherInfoSize = 0;
53 }
54
getTile(const ossimIrect & tileRect,ossim_uint32 resLevel)55 ossimRefPtr<ossimImageData> ossimFeatherMosaic::getTile(const ossimIrect& tileRect,
56 ossim_uint32 resLevel)
57 {
58 long w = tileRect.width();
59 long h = tileRect.height();
60 ossimIpt origin = tileRect.ul();
61
62 if(!isSourceEnabled())
63 {
64 return ossimImageMosaic::getTile(tileRect, resLevel);
65 }
66 if(!theTile||!theAlphaSum||!theResult||!theInputFeatherInformation)
67 {
68 initialize();
69
70 if(!theTile||!theAlphaSum||!theResult||!theInputFeatherInformation)
71 {
72 return ossimImageMosaic::getTile(tileRect, resLevel);
73 }
74 }
75 ossim_uint32 size = getNumberOfInputs();
76 theAlphaSum->setImageRectangle(tileRect);
77 theResult->setImageRectangle(tileRect);
78
79 if(size == 0)
80 {
81 return ossimRefPtr<ossimImageData>();
82 }
83
84 if(size == 1)
85 {
86 return ossimImageMosaic::getTile(tileRect, resLevel);
87 }
88
89 long tileW = theTile->getWidth();
90 long tileH = theTile->getHeight();
91 if((w != tileW)||
92 (h != tileH))
93 {
94 theTile->setWidth(w);
95 theTile->setHeight(h);
96 if((w*h)!=(tileW*tileH))
97 {
98 theTile->initialize();
99 }
100 }
101 theTile->setOrigin(origin);
102 theTile->makeBlank();
103
104 switch(theTile->getScalarType())
105 {
106 case OSSIM_UCHAR:
107 {
108 return combine(static_cast<ossim_uint8>(0),
109 tileRect, resLevel);
110 }
111 case OSSIM_USHORT16:
112 case OSSIM_USHORT11:
113 case OSSIM_USHORT12:
114 case OSSIM_USHORT13:
115 case OSSIM_USHORT14:
116 case OSSIM_USHORT15:
117 {
118 return combine(static_cast<ossim_uint16>(0),
119 tileRect, resLevel);
120 }
121 case OSSIM_SSHORT16:
122 {
123 return combine(static_cast<ossim_sint16>(0),
124 tileRect, resLevel);
125 }
126 case OSSIM_DOUBLE:
127 case OSSIM_NORMALIZED_DOUBLE:
128 {
129 return combine(static_cast<double>(0),
130 tileRect, resLevel);
131 }
132 case OSSIM_FLOAT:
133 case OSSIM_NORMALIZED_FLOAT:
134 {
135 return combine(static_cast<float>(0),
136 tileRect, resLevel);
137 }
138 case OSSIM_SCALAR_UNKNOWN:
139 default:
140 {
141 ossimNotify(ossimNotifyLevel_WARN)
142 << "ossimFeatherMosaic::getTile: error, unknown scalar type!!!"
143 << std::endl;
144 }
145 }
146 return ossimRefPtr<ossimImageData>();
147 }
148
149
combine(T,const ossimIrect & tileRect,ossim_uint32 resLevel)150 template <class T> ossimRefPtr<ossimImageData> ossimFeatherMosaic::combine(
151 T,
152 const ossimIrect& tileRect,
153 ossim_uint32 resLevel)
154 {
155 ossimRefPtr<ossimImageData> currentImageData;
156 ossim_uint32 band;
157 long upperBound = theTile->getWidth()*theTile->getHeight();
158 long offset = 0;
159 long row = 0;
160 long col = 0;
161 long numberOfTilesProcessed = 0;
162 float *sumBand = static_cast<float*>(theAlphaSum->getBuf());
163 float *bandRes = NULL;
164 ossimIpt point;
165
166 theAlphaSum->fill(0.0);
167 theResult->fill(0.0);
168
169 T** srcBands = new T*[theLargestNumberOfInputBands];
170 T** destBands = new T*[theLargestNumberOfInputBands];
171 ossim_uint32 layerIdx = 0;
172 currentImageData = getNextTile(layerIdx,
173 0,
174 tileRect,
175 resLevel);
176 if(!currentImageData.valid())
177 {
178 delete [] srcBands;
179 delete [] destBands;
180 return currentImageData;
181 }
182 ossim_uint32 minNumberOfBands = currentImageData->getNumberOfBands();
183 for(band = 0; band < minNumberOfBands; ++band)
184 {
185 srcBands[band] = static_cast<T*>(currentImageData->getBuf(band));
186 destBands[band] = static_cast<T*>(theTile->getBuf(band));
187 }
188 // if the src is smaller than the destination in number
189 // of bands we will just duplicate the last band.
190 for(;band < theLargestNumberOfInputBands; ++band)
191 {
192 srcBands[band] = srcBands[minNumberOfBands - 1];
193 destBands[band] = static_cast<T*>(theTile->getBuf(band));
194 }
195
196 // most of the time we will not overlap so let's
197 // copy the first tile into destination and check later.
198 //
199 for(band = 0; band < theTile->getNumberOfBands();++band)
200 {
201 T* destBand = destBands[band];
202 T* srcBand = srcBands[band];
203 if(destBand&&srcBand)
204 {
205 for(offset = 0; offset < upperBound;++offset)
206 {
207 *destBand = *srcBand;
208 ++srcBand; ++destBand;
209 }
210 }
211 }
212 theTile->setDataObjectStatus(currentImageData->getDataObjectStatus());
213
214 while(currentImageData.valid())
215 {
216 ossimDataObjectStatus currentStatus = currentImageData->getDataObjectStatus();
217 point = currentImageData->getOrigin();
218 long h = (long)currentImageData->getHeight();
219 long w = (long)currentImageData->getWidth();
220 if( (currentStatus != OSSIM_EMPTY) &&
221 (currentStatus != OSSIM_NULL))
222 {
223 ++numberOfTilesProcessed;
224 offset = 0;
225 minNumberOfBands = currentImageData->getNumberOfBands();
226 for(band = 0; band < minNumberOfBands; ++band)
227 {
228 srcBands[band] = static_cast<T*>(currentImageData->getBuf(band));
229 }
230 // if the src is smaller than the destination in number
231 // of bands we will just duplicate the last band.
232 for(;band < theLargestNumberOfInputBands; ++band)
233 {
234 srcBands[band] = srcBands[minNumberOfBands - 1];
235 }
236 if(currentStatus == OSSIM_PARTIAL)
237 {
238 for(row = 0; row < h; ++row)
239 {
240 for(col = 0; col < w; ++col)
241 {
242 if(!currentImageData->isNull(offset))
243 {
244 double weight = computeWeight(layerIdx,
245 ossimDpt(point.x+col,
246 point.y+row));
247
248 for(band = 0; band < theLargestNumberOfInputBands; ++band)
249 {
250 bandRes = static_cast<float*>(theResult->getBuf(band));
251 bandRes[offset] += (srcBands[band][offset]*weight);
252 }
253 sumBand[offset] += weight;
254 }
255 ++offset;
256 }
257 }
258 }
259 else
260 {
261 offset = 0;
262
263 for(row = 0; row < h; ++row)
264 {
265 for(col = 0; col < w; ++col)
266 {
267 double weight = computeWeight(layerIdx,
268 ossimDpt(point.x+col,
269 point.y+row));
270
271 for(band = 0; band < theLargestNumberOfInputBands; ++band)
272 {
273 bandRes = static_cast<float*>(theResult->getBuf(band));
274
275 bandRes[offset] += (srcBands[band][offset]*weight);
276 }
277 sumBand[offset] += weight;
278 ++offset;
279 }
280 }
281 }
282 }
283 currentImageData = getNextTile(layerIdx, tileRect, resLevel);
284 }
285 upperBound = theTile->getWidth()*theTile->getHeight();
286
287 if(numberOfTilesProcessed > 1)
288 {
289 const double* minPix = theTile->getMinPix();
290 const double* maxPix = theTile->getMaxPix();
291 const double* nullPix= theTile->getNullPix();
292 for(offset = 0; offset < upperBound;++offset)
293 {
294 for(band = 0; band < theTile->getNumberOfBands();++band)
295 {
296 T* destBand = static_cast<T*>(theTile->getBuf(band));
297 float* weightedBand = static_cast<float*>(theResult->getBuf(band));
298
299 // this should be ok to test 0.0 instead of
300 // FLT_EPSILON range for 0 since we set it.
301 if(sumBand[offset] != 0.0)
302 {
303 weightedBand[offset] = (weightedBand[offset])/sumBand[offset];
304 if(weightedBand[offset]<minPix[band])
305 {
306 weightedBand[offset] = minPix[band];
307 }
308 else if(weightedBand[offset] > maxPix[band])
309 {
310 weightedBand[offset] = maxPix[band];
311 }
312 }
313 else
314 {
315 weightedBand[offset] = nullPix[band];
316 }
317 destBand[offset] = static_cast<T>(weightedBand[offset]);
318 }
319 }
320 theTile->validate();
321 }
322
323 delete [] srcBands;
324 delete [] destBands;
325
326 return theTile;
327 }
328
computeWeight(long index,const ossimDpt & point) const329 double ossimFeatherMosaic::computeWeight(long index,
330 const ossimDpt& point)const
331 {
332 ossimFeatherInputInformation& info = theInputFeatherInformation[index];
333 double result = 0.0;
334 ossimDpt delta = point-info.theCenter;
335
336 double length1 = fabs(delta.x*info.theAxis1.x + delta.y*info.theAxis1.y)/info.theAxis1Length;
337 double length2 = fabs(delta.x*info.theAxis2.x + delta.y*info.theAxis2.y)/info.theAxis2Length;
338
339 if(length1 > length2)
340 {
341 result = (1.0 - length1);
342 }
343 else
344 {
345 result = (1.0 - length2);
346 }
347 if(result < 0) result = 0;
348
349 return result;
350 }
351
initialize()352 void ossimFeatherMosaic::initialize()
353 {
354 ossimImageMosaic::initialize();
355
356 allocate();
357 if(theTile.valid())
358 {
359 theAlphaSum = new ossimImageData(this,
360 OSSIM_FLOAT,
361 1,
362 theTile->getWidth(),
363 theTile->getHeight());
364 theResult = new ossimImageData(this,
365 OSSIM_FLOAT,
366 theLargestNumberOfInputBands,
367 theTile->getWidth(),
368 theTile->getHeight());
369 theAlphaSum->initialize();
370 theResult->initialize();
371 }
372 std::vector<ossimIpt> validVertices;
373 if(!getNumberOfInputs()) return;
374
375 if(theInputFeatherInformation)
376 {
377 delete [] theInputFeatherInformation;
378 theInputFeatherInformation = NULL;
379 }
380 theFeatherInfoSize = getNumberOfInputs();
381 theInputFeatherInformation = new ossimFeatherInputInformation[theFeatherInfoSize];
382 for(long index = 0; index < theFeatherInfoSize; ++ index)
383 {
384 validVertices.clear();
385 ossimImageSource* temp = PTR_CAST(ossimImageSource, getInput(index));
386 if(temp)
387 {
388 temp->getValidImageVertices(validVertices, OSSIM_CLOCKWISE_ORDER);
389 theInputFeatherInformation[index].setVertexList(validVertices);
390 }
391 }
392 }
393
394
ossimFeatherInputInformation(const std::vector<ossimIpt> & validVertices)395 ossimFeatherMosaic::ossimFeatherInputInformation::ossimFeatherInputInformation(const std::vector<ossimIpt>& validVertices)
396 {
397 setVertexList(validVertices);
398 }
399
400
setVertexList(const std::vector<ossimIpt> & validVertices)401 void ossimFeatherMosaic::ossimFeatherInputInformation::setVertexList(const std::vector<ossimIpt>& validVertices)
402 {
403 const char* MODULE = "ossimFeatherMosaic::ossimFeatherInputInformation::setVertexList()";
404
405 theValidVertices = validVertices;
406
407 theCenter = ossimDpt(0,0);
408 theAxis1 = ossimDpt(1, 0);
409 theAxis2 = ossimDpt(0, 1);
410 theAxis1Length = 1;
411 theAxis2Length = 1;
412
413 double xSum=0.0, ySum=0.0;
414 ossim_uint32 upperBound = (ossim_uint32)validVertices.size();
415 if(upperBound)
416 {
417 for(ossim_uint32 index = 0; index < upperBound; ++index)
418 {
419 xSum += validVertices[index].x;
420 ySum += validVertices[index].y;
421 }
422
423 theCenter.x = xSum/upperBound;
424 theCenter.y = ySum/upperBound;
425
426 // for now we just want a quick implementation of something
427 // and we know that we have 4 vertices for the bounding valid
428 // vertices.
429 //
430 if(upperBound == 4)
431 {
432 ossimDpt edgeDirection1 = validVertices[1] - validVertices[0];
433 ossimDpt edgeDirection2 = validVertices[2] - validVertices[1];
434
435 theAxis1 = ossimDpt(-edgeDirection1.y, edgeDirection1.x);
436
437 theAxis2 = ossimDpt(-edgeDirection2.y, edgeDirection2.x);
438
439 theAxis1 = theAxis1/theAxis1.length();
440 theAxis2 = theAxis2/theAxis2.length();
441
442 ossimLine line1(theCenter,
443 theCenter + theAxis1*2);
444 ossimLine line2(validVertices[1],
445 validVertices[0]);
446 ossimLine line3(theCenter,
447 theCenter + theAxis2*2);
448 ossimLine line4(validVertices[2],
449 validVertices[1]);
450
451 ossimDpt intersectionPoint1 = line1.intersectInfinite(line2);
452 ossimDpt intersectionPoint2 = line3.intersectInfinite(line4);
453
454
455 theAxis1Length = ossim::round<int>((theCenter-intersectionPoint1).length());
456 theAxis2Length = ossim::round<int>((theCenter-intersectionPoint2).length());
457
458 if(traceDebug())
459 {
460 CLOG << "theAxis1Length: " << theAxis1Length << endl
461 << "theAxis2Length: " << theAxis2Length << endl
462 << "center: " << theCenter << endl;
463 }
464 }
465 }
466 }
467
operator <<(ostream & out,const ossimFeatherMosaic::ossimFeatherInputInformation & data)468 ostream& operator<<(ostream& out,
469 const ossimFeatherMosaic::ossimFeatherInputInformation& data)
470 {
471 out << "center: " << data.theCenter << endl
472 << "axis1: " << data.theAxis1 << endl
473 << "axis2: " << data.theAxis2 << endl
474 << "axis1_length: " << data.theAxis1Length << endl
475 << "axis2_length: " << data.theAxis2Length << endl
476 << "valid vertices: " << endl;
477 std::copy(data.theValidVertices.begin(),
478 data.theValidVertices.end(),
479 std::ostream_iterator<ossimDpt>(out, "\n"));
480 return out;
481 }
482